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