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