ifa: Doc some heuristics about EAP quadrature
[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 /**
892  * MAP is not affected by the number of items. EAP is. Likelihood can
893  * get concentrated in a single quadrature ordinate. For 3PL, response
894  * patterns can have a bimodal likelihood. This will confuse MAP and
895  * is a key advantage of EAP (Thissen & Orlando, 2001, p. 136).
896  *
897  * Thissen, D. & Orlando, M. (2001). IRT for items scored in two
898  * categories. In D. Thissen & H. Wainer (Eds.), \emph{Test scoring}
899  * (pp 73-140). Lawrence Erlbaum Associates, Inc.
900  */
901 omxRListElement *
902 ba81EAP(omxExpectation *oo, int *numReturns)
903 {
904         omxBA81State *state = (omxBA81State *) oo->argStruct;
905         int maxDims = state->maxDims;
906         //int numSpecific = state->numSpecific;
907
908         *numReturns = 2; // + (maxDims > 1) + (numSpecific > 1);
909         omxRListElement *out = (omxRListElement*) R_alloc(*numReturns, sizeof(omxRListElement));
910
911         out[0].numValues = 1;
912         out[0].values = (double*) R_alloc(1, sizeof(double));
913         strcpy(out[0].label, "Minus2LogLikelihood");
914         out[0].values[0] = state->ll;
915
916         omxData *data = state->data;
917         int numUnique = state->numUnique;
918
919         // TODO Wainer & Thissen. (1987). Estimating ability with the wrong
920         // model. Journal of Educational Statistics, 12, 339-368.
921
922         int numQpoints = state->numQpoints * 2;  // make configurable TODO
923
924         if (numQpoints < 1 + 2.0 * sqrt(state->itemSpec->cols)) {
925                 // Thissen & Orlando (2001, p. 136)
926                 warning("EAP requires at least 2*sqrt(items) quadrature points");
927         }
928
929         double Qpoint[numQpoints];
930         double Qarea[numQpoints];
931         const double Qwidth = 4;
932         for (int qx=0; qx < numQpoints; qx++) {
933                 Qpoint[qx] = Qwidth - qx * Qwidth*2 / (numQpoints-1);
934                 Qarea[qx] = 1.0/numQpoints;
935         }
936         ba81SetupQuadrature(oo, numQpoints, Qpoint, Qarea);
937         ba81Estep(oo);   // recalc patternLik with a flat prior
938
939         double *cov = NULL;
940         /*
941         if (maxDims > 1) {
942                 strcpy(out[2].label, "ability.cov");
943                 out[2].numValues = -1;
944                 out[2].rows = maxDims;
945                 out[2].cols = maxDims;
946                 out[2].values = (double*) R_alloc(out[2].rows * out[2].cols, sizeof(double));
947                 cov = out[2].values;
948                 OMXZERO(cov, out[2].rows * out[2].cols);
949         }
950         */
951         double *spstats = NULL;
952         /*
953         if (numSpecific) {
954                 strcpy(out[3].label, "specific");
955                 out[3].numValues = -1;
956                 out[3].rows = numSpecific;
957                 out[3].cols = 2;
958                 out[3].values = (double*) R_alloc(out[3].rows * out[3].cols, sizeof(double));
959                 spstats = out[3].values;
960         }
961         */
962
963         // allocation of workspace could be optional
964         int numThreads = getNumThreads(oo);
965         double *workspace = Realloc(NULL, numUnique * maxDims * 2 * numThreads, double);
966
967         // Need a separate work space because the destination needs
968         // to be in unsorted order with duplicated rows.
969         double *ability = Calloc(numUnique * maxDims * 2, double);
970
971 #pragma omp parallel for num_threads(oo->currentState->numThreads)
972         for (long qx=0; qx < state->totalQuadPoints; qx++) {
973                 ba81EAP1(oo, workspace, qx, maxDims, numUnique, ability, cov, spstats);
974         }
975
976         /*
977         // make symmetric
978         for (int d1=0; d1 < maxDims; d1++) {
979                 for (int d2=0; d2 < d1; d2++) {
980                         cov[d2 * maxDims + d1] = cov[d1 * maxDims + d2];
981                 }
982         }
983         */
984
985 #pragma omp parallel for num_threads(oo->currentState->numThreads)
986         for (long qx=0; qx < state->totalQuadPoints; qx++) {
987                 ba81EAP2(oo, workspace, qx, maxDims, numUnique, ability, spstats);
988         }
989
990         for (int px=0; px < numUnique; px++) {
991                 double *arow = ability + px * 2 * maxDims;
992                 for (int dx=0; dx < maxDims; dx++) {
993                         arow[dx*2+1] = sqrt(arow[dx*2+1]);
994                 }
995         }
996
997         strcpy(out[1].label, "ability");
998         out[1].numValues = -1;
999         out[1].rows = data->rows;
1000         out[1].cols = 2 * maxDims;
1001         out[1].values = (double*) R_alloc(out[1].rows * out[1].cols, sizeof(double));
1002
1003         for (int rx=0; rx < numUnique; rx++) {
1004                 double *pa = ability + rx * 2 * maxDims;
1005
1006                 int dups = omxDataNumIdenticalRows(state->data, state->rowMap[rx]);
1007                 for (int dup=0; dup < dups; dup++) {
1008                         int dest = omxDataIndex(data, state->rowMap[rx]+dup);
1009                         int col=0;
1010                         for (int dx=0; dx < maxDims; dx++) {
1011                                 out[1].values[col * out[1].rows + dest] = pa[col]; ++col;
1012                                 out[1].values[col * out[1].rows + dest] = pa[col]; ++col;
1013                         }
1014                 }
1015         }
1016         Free(ability);
1017         Free(workspace);
1018         return out;
1019 }
1020
1021 static void ba81Destroy(omxExpectation *oo) {
1022         if(OMX_DEBUG) {
1023                 Rprintf("Freeing %s function.\n", NAME);
1024         }
1025         omxBA81State *state = (omxBA81State *) oo->argStruct;
1026         Rprintf("fit %d gradient %d\n", state->fitCount, state->gradientCount);
1027         omxFreeAllMatrixData(state->itemSpec);
1028         omxFreeAllMatrixData(state->itemParam);
1029         omxFreeAllMatrixData(state->EitemParam);
1030         omxFreeAllMatrixData(state->design);
1031         omxFreeAllMatrixData(state->customPrior);
1032         Free(state->logNumIdentical);
1033         Free(state->Qpoint);
1034         Free(state->Qarea);
1035         Free(state->logQarea);
1036         Free(state->rowMap);
1037         Free(state->patternLik);
1038         Free(state->lxk);
1039         Free(state->Slxk);
1040         Free(state->allSlxk);
1041         Free(state->Sgroup);
1042         Free(state->paramMap);
1043         Free(state);
1044 }
1045
1046 int ba81ExpectationHasGradients(omxExpectation* oo)
1047 {
1048         omxBA81State *state = (omxBA81State *) oo->argStruct;
1049         return state->computeRPF == standardComputeRPF;
1050 }
1051
1052 void omxInitExpectationBA81(omxExpectation* oo) {
1053         omxState* currentState = oo->currentState;      
1054         SEXP rObj = oo->rObj;
1055         SEXP tmp;
1056         
1057         if(OMX_DEBUG) {
1058                 Rprintf("Initializing %s.\n", NAME);
1059         }
1060         if (!rpf_model) {
1061                 get_librpf_t get_librpf = (get_librpf_t) R_GetCCallable("rpf", "get_librpf_model");
1062                 (*get_librpf)(&rpf_numModels, &rpf_model);
1063         }
1064         
1065         omxBA81State *state = Calloc(1, omxBA81State);
1066         oo->argStruct = (void*) state;
1067
1068         state->ll = 1e20;   // finite but big
1069         
1070         PROTECT(tmp = GET_SLOT(rObj, install("data")));
1071         state->data = omxNewDataFromMxDataPtr(tmp, currentState);
1072         UNPROTECT(1);
1073
1074         if (strcmp(omxDataType(state->data), "raw") != 0) {
1075                 omxRaiseErrorf(currentState, "%s unable to handle data type %s", NAME, omxDataType(state->data));
1076                 return;
1077         }
1078
1079         PROTECT(state->rpf = GET_SLOT(rObj, install("RPF")));
1080         if (state->rpf == R_NilValue) {
1081                 state->computeRPF = standardComputeRPF;
1082         } else {
1083                 state->computeRPF = RComputeRPF;
1084         }
1085
1086         state->itemSpec =
1087                 omxNewMatrixFromIndexSlot(rObj, currentState, "ItemSpec");
1088         state->design =
1089                 omxNewMatrixFromIndexSlot(rObj, currentState, "Design");
1090         state->itemParam =
1091                 omxNewMatrixFromIndexSlot(rObj, currentState, "ItemParam");
1092         state->EitemParam =
1093                 omxInitTemporaryMatrix(NULL, state->itemParam->rows, state->itemParam->cols,
1094                                        TRUE, currentState);
1095         state->customPrior =
1096                 omxNewMatrixFromIndexSlot(rObj, currentState, "CustomPrior");
1097         
1098         oo->computeFun = ba81Estep;
1099         oo->destructFun = ba81Destroy;
1100         
1101         // TODO: Exactly identical rows do not contribute any information.
1102         // The sorting algorithm ought to remove them so we don't waste RAM.
1103         // The following summary stats would be cheaper to calculate too.
1104
1105         int numUnique = 0;
1106         omxData *data = state->data;
1107         if (omxDataNumFactor(data) != data->cols) {
1108                 // verify they are ordered factors TODO
1109                 omxRaiseErrorf(currentState, "%s: all columns must be factors", NAME);
1110                 return;
1111         }
1112
1113         for (int rx=0; rx < data->rows;) {
1114                 rx += omxDataNumIdenticalRows(state->data, rx);
1115                 ++numUnique;
1116         }
1117         state->numUnique = numUnique;
1118
1119         state->rowMap = Realloc(NULL, numUnique, int);
1120         state->logNumIdentical = Realloc(NULL, numUnique, double);
1121
1122         int numItems = state->itemParam->cols;
1123
1124         for (int rx=0, ux=0; rx < data->rows; ux++) {
1125                 if (rx == 0) {
1126                         // all NA rows will sort to the top
1127                         int na=0;
1128                         for (int ix=0; ix < numItems; ix++) {
1129                                 if (omxIntDataElement(data, 0, ix) == NA_INTEGER) { ++na; }
1130                         }
1131                         if (na == numItems) {
1132                                 omxRaiseErrorf(currentState, "Remove rows with all NAs");
1133                                 return;
1134                         }
1135                 }
1136                 int dups = omxDataNumIdenticalRows(state->data, rx);
1137                 state->logNumIdentical[ux] = log(dups);
1138                 state->rowMap[ux] = rx;
1139                 rx += dups;
1140         }
1141
1142         state->patternLik = Realloc(NULL, numUnique, double);
1143
1144         int numThreads = getNumThreads(oo);
1145
1146         int maxSpec = 0;
1147         int maxParam = 0;
1148         state->maxDims = 0;
1149         state->maxOutcomes = 0;
1150
1151         for (int cx = 0; cx < data->cols; cx++) {
1152                 const double *spec = omxMatrixColumn(state->itemSpec, cx);
1153                 int id = spec[RPF_ISpecID];
1154                 if (id < 0 || id >= rpf_numModels) {
1155                         omxRaiseErrorf(currentState, "ItemSpec column %d has unknown item model %d", cx, id);
1156                         return;
1157                 }
1158
1159                 int dims = spec[RPF_ISpecDims];
1160                 if (state->maxDims < dims)
1161                         state->maxDims = dims;
1162
1163                 // TODO verify that item model can have requested number of outcomes
1164                 int no = spec[RPF_ISpecOutcomes];
1165                 if (state->maxOutcomes < no)
1166                         state->maxOutcomes = no;
1167
1168                 int numSpec = (*rpf_model[id].numSpec)(spec);
1169                 if (maxSpec < numSpec)
1170                         maxSpec = numSpec;
1171
1172                 int numParam = (*rpf_model[id].numParam)(spec);
1173                 if (maxParam < numParam)
1174                         maxParam = numParam;
1175         }
1176
1177         if (state->itemSpec->cols != data->cols || state->itemSpec->rows != maxSpec) {
1178                 omxRaiseErrorf(currentState, "ItemSpec must have %d item columns and %d rows",
1179                                data->cols, maxSpec);
1180                 return;
1181         }
1182         if (state->itemParam->rows != maxParam) {
1183                 omxRaiseErrorf(currentState, "ItemParam should have %d rows", maxParam);
1184                 return;
1185         }
1186
1187         if (state->design == NULL) {
1188                 state->maxAbilities = state->maxDims;
1189                 state->design = omxInitTemporaryMatrix(NULL, state->maxDims, numItems,
1190                                        TRUE, currentState);
1191                 for (int ix=0; ix < numItems; ix++) {
1192                         for (int dx=0; dx < state->maxDims; dx++) {
1193                                 omxSetMatrixElement(state->design, dx, ix, (double)dx+1);
1194                         }
1195                 }
1196         } else {
1197                 omxMatrix *design = state->design;
1198                 if (design->cols != numItems ||
1199                     design->rows != state->maxDims) {
1200                         omxRaiseErrorf(currentState, "Design matrix should have %d rows and %d columns",
1201                                        state->maxDims, numItems);
1202                         return;
1203                 }
1204
1205                 state->maxAbilities = 0;
1206                 for (int ix=0; ix < design->rows * design->cols; ix++) {
1207                         double got = design->data[ix];
1208                         if (!R_FINITE(got)) continue;
1209                         if (round(got) != got) error("Design matrix can only contain integers"); // TODO better way?
1210                         if (state->maxAbilities < got)
1211                                 state->maxAbilities = got;
1212                 }
1213         }
1214         if (state->maxAbilities <= state->maxDims) {
1215                 state->Sgroup = Calloc(numItems, int);
1216         } else {
1217                 int Sgroup0 = state->maxDims;
1218                 state->Sgroup = Realloc(NULL, numItems, int);
1219                 for (int ix=0; ix < numItems; ix++) {
1220                         int ss=-1;
1221                         for (int dx=0; dx < state->maxDims; dx++) {
1222                                 int ability = omxMatrixElement(state->design, dx, ix);
1223                                 if (ability >= Sgroup0) {
1224                                         if (ss == -1) {
1225                                                 ss = ability;
1226                                         } else {
1227                                                 omxRaiseErrorf(currentState, "Item %d cannot belong to more than "
1228                                                                "1 specific dimension (both %d and %d)",
1229                                                                ix, ss, ability);
1230                                                 return;
1231                                         }
1232                                 }
1233                         }
1234                         if (ss == -1) ss = 0;
1235                         state->Sgroup[ix] = ss - Sgroup0;
1236                 }
1237                 state->numSpecific = state->maxAbilities - state->maxDims + 1;
1238                 state->allSlxk = Realloc(NULL, numUnique * numThreads, double);
1239                 state->Slxk = Realloc(NULL, numUnique * state->numSpecific * numThreads, double);
1240         }
1241
1242         PROTECT(tmp = GET_SLOT(rObj, install("cache")));
1243         state->cacheLXK = asLogical(tmp);
1244
1245         PROTECT(tmp = GET_SLOT(rObj, install("GHpoints")));
1246         double *qpoints = REAL(tmp);
1247         int numQPoints = length(tmp);
1248
1249         PROTECT(tmp = GET_SLOT(rObj, install("GHarea")));
1250         double *qarea = REAL(tmp);
1251         if (numQPoints != length(tmp)) error("length(GHpoints) != length(GHarea)");
1252
1253         ba81SetupQuadrature(oo, numQPoints, qpoints, qarea);
1254
1255         // verify data bounded between 1 and numOutcomes TODO
1256         // hm, looks like something could be added to omxData for column summary stats?
1257 }