Revive latent distribution gradients
[openmx:openmx.git] / src / omxFitFunctionBA81.cpp
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 #include <algorithm>
19
20 #include "omxFitFunction.h"
21 #include "omxExpectationBA81.h"
22 #include "omxOpenmpWrap.h"
23 #include "libifa-rpf.h"
24 #include "matrix.h"
25 #include "omxBuffer.h"
26
27 struct BA81FitState {
28
29         int haveLatentMap;
30         std::vector<int> latentMap;
31         bool freeLatents;
32         int ElatentVersion;
33
34         int haveItemMap;
35         size_t numFreeParam;
36         int itemDerivPadSize;     // maxParam + maxParam*(1+maxParam)/2
37         std::vector<int> paramPerItem;       // itemParam->cols
38         std::vector<int> paramFlavor;        // freeParam
39         std::vector<int> paramMap;           // itemParam->cols * itemDerivPadSize -> index of free parameter
40         std::vector<int> itemGradMap;        // index of gradient -> index of free parameter
41         std::vector<int> paramLocations;     // param# -> count of appearances in ItemParam
42         std::vector<int> itemParamFree;      // itemParam->cols * itemParam->rows
43         std::vector<int> ihessDivisor;       // freeParam * freeParam
44         std::vector< matrixVectorProdTerm > hgProd;
45
46         omxMatrix *itemParam;
47         omxMatrix *latentMean;
48         omxMatrix *latentCov;
49
50         BA81FitState();
51         ~BA81FitState();
52         void copyEstimates(BA81Expect *estate);
53 };
54
55 // writes to upper triangle of full matrix
56 static void addSymOuterProd(const double weight, const double *vec, const int len, double *out)
57 {
58         for (int d1=0; d1 < len; ++d1) {
59                 for (int d2=0; d2 <= d1; ++d2) {
60                         out[d1 * len + d2] += weight * vec[d1] * vec[d2];
61                 }
62         }
63 }
64
65 BA81FitState::BA81FitState()
66 {
67         haveItemMap = FREEVARGROUP_INVALID;
68         haveLatentMap = FREEVARGROUP_INVALID;
69 }
70
71 void BA81FitState::copyEstimates(BA81Expect *estate)
72 {
73         omxCopyMatrix(itemParam, estate->itemParam);
74         omxCopyMatrix(latentMean, estate->latentMeanOut);
75         omxCopyMatrix(latentCov, estate->latentCovOut);
76 }
77
78 static void buildLatentParamMap(omxFitFunction* oo, FitContext *fc)
79 {
80         FreeVarGroup *fvg = fc->varGroup;
81         BA81FitState *state = (BA81FitState *) oo->argStruct;
82         std::vector<int> &latentMap = state->latentMap;
83         BA81Expect *estate = (BA81Expect*) oo->expectation->argStruct;
84         int meanNum = estate->latentMeanOut->matrixNumber;
85         int covNum = estate->latentCovOut->matrixNumber;
86         int maxAbilities = estate->maxAbilities;
87         int numLatents = maxAbilities + triangleLoc1(maxAbilities);
88
89         if (state->haveLatentMap == fc->varGroup->id) return;
90         if (estate->verbose) mxLog("%s: rebuild latent parameter map for %d", oo->matrix->name, fc->varGroup->id);
91
92         state->freeLatents = false;
93         latentMap.assign(numLatents, -1);
94
95         int numParam = int(fvg->vars.size());
96         for (int px=0; px < numParam; px++) {
97                 omxFreeVar *fv = fvg->vars[px];
98                 for (size_t lx=0; lx < fv->locations.size(); lx++) {
99                         omxFreeVarLocation *loc = &fv->locations[lx];
100                         int matNum = ~loc->matrix;
101                         if (matNum == meanNum) {
102                                 latentMap[loc->row + loc->col] = px;
103                                 state->freeLatents = true;
104                         } else if (matNum == covNum) {
105                                 int a1 = loc->row;
106                                 int a2 = loc->col;
107                                 if (a1 < a2) std::swap(a1, a2);
108                                 int cell = maxAbilities + triangleLoc1(a1) + a2;
109                                 if (latentMap[cell] == -1) {
110                                         latentMap[cell] = px;
111
112                                         if (a1 == a2 && fv->lbound == NEG_INF) {
113                                                 fv->lbound = 1e-6;  // variance must be positive
114                                                 if (fc->est[px] < fv->lbound) {
115                                                         error("Starting value for variance %s is negative", fv->name);
116                                                 }
117                                         }
118                                 } else if (latentMap[cell] != px) {
119                                         // doesn't detect similar problems in multigroup constraints TODO
120                                         error("In covariance matrix, %s and %s must be constrained equal to preserve symmetry",
121                                               fvg->vars[latentMap[cell]]->name, fv->name);
122                                 }
123                                 state->freeLatents = true;
124                         }
125                 }
126         }
127         state->haveLatentMap = fc->varGroup->id;
128 }
129
130 static void buildItemParamMap(omxFitFunction* oo, FitContext *fc)
131 {
132         FreeVarGroup *fvg = fc->varGroup;
133         BA81FitState *state = (BA81FitState *) oo->argStruct;
134         BA81Expect *estate = (BA81Expect*) oo->expectation->argStruct;
135
136         if (state->haveItemMap == fc->varGroup->id) return;
137         if (estate->verbose) mxLog("%s: rebuild item parameter map for %d", oo->matrix->name, fc->varGroup->id);
138
139         omxMatrix *itemParam = estate->itemParam;
140         int size = itemParam->cols * state->itemDerivPadSize;
141         state->paramMap.assign(size, -1);  // matrix location to free param index
142         state->itemParamFree.assign(itemParam->rows * itemParam->cols, FALSE);
143
144         size_t numFreeParams = state->numFreeParam = fvg->vars.size();
145         state->paramLocations.assign(numFreeParams, 0);
146         state->paramFlavor.assign(numFreeParams, -1);
147
148         int totalParam = 0;
149         state->paramPerItem.resize(itemParam->cols);
150         for (int cx=0; cx < itemParam->cols; ++cx) {
151                 const double *spec = estate->itemSpec[cx];
152                 const int id = spec[RPF_ISpecID];
153                 const int numParam = (*rpf_model[id].numParam)(spec);
154                 state->paramPerItem[cx] = numParam;
155                 totalParam += numParam;
156         }
157         state->itemGradMap.assign(totalParam, -1);
158
159         for (size_t px=0; px < numFreeParams; px++) {
160                 omxFreeVar *fv = fvg->vars[px];
161                 state->paramLocations[px] = int(fv->locations.size());
162                 for (size_t lx=0; lx < fv->locations.size(); lx++) {
163                         omxFreeVarLocation *loc = &fv->locations[lx];
164                         int matNum = ~loc->matrix;
165                         if (matNum == itemParam->matrixNumber) {
166                                 int at = loc->col * state->itemDerivPadSize + loc->row;
167                                 state->paramMap[at] = px;
168                                 state->itemParamFree[loc->col * itemParam->rows + loc->row] = TRUE;
169
170                                 const double *spec = estate->itemSpec[loc->col];
171                                 int id = spec[RPF_ISpecID];
172                                 int flavor;
173                                 double upper, lower;
174                                 (*rpf_model[id].paramInfo)(spec, loc->row, &flavor, &upper, &lower);
175                                 if (state->paramFlavor[px] < 0) {
176                                         state->paramFlavor[px] = flavor;
177                                 } else if (state->paramFlavor[px] != flavor) {
178                                         error("Cannot equate %s with %s[%d,%d]", fv->name,
179                                               itemParam->name, loc->row, loc->col);
180                                 }
181                                 if (fv->lbound == NEG_INF && isfinite(lower)) {
182                                         fv->lbound = lower;
183                                         if (fc->est[px] < fv->lbound) {
184                                                 error("Starting value %s %f less than lower bound %f",
185                                                       fv->name, fc->est[px], lower);
186                                         }
187                                 }
188                                 if (fv->ubound == INF && isfinite(upper)) {
189                                         fv->ubound = upper;
190                                         if (fc->est[px] > fv->ubound) {
191                                                 error("Starting value %s %f greater than upper bound %f",
192                                                       fv->name, fc->est[px], upper);
193                                         }
194                                 }
195                         }
196                 }
197         }
198
199         int gradOffset = 0;
200         for (int cx=0; cx < itemParam->cols; ++cx) {
201                 for (int rx=0; rx < state->paramPerItem[cx]; ++rx) {
202                         int at = cx * state->itemDerivPadSize + rx;
203                         int px = state->paramMap[at];
204                         if (px >= 0) state->itemGradMap[gradOffset] = px;
205                         ++gradOffset;
206                 }
207         }
208
209         state->ihessDivisor.resize(size);
210
211         for (int cx=0; cx < itemParam->cols; ++cx) {
212                 int numParam = state->paramPerItem[cx];
213                 for (int p1=0; p1 < numParam; p1++) {
214                         int at1 = state->paramMap[cx * state->itemDerivPadSize + p1];
215                         if (at1 < 0) continue;
216
217                         for (int p2=0; p2 <= p1; p2++) {
218                                 int at2 = state->paramMap[cx * state->itemDerivPadSize + p2];
219                                 if (at2 < 0) continue;
220
221                                 if (at1 < at2) std::swap(at1, at2);  // lower triangle
222
223                                 //mxLog("Item %d param(%d,%d) -> H[%d,%d]", cx, p1, p2, at1, at2);
224                                 int at = cx * state->itemDerivPadSize + numParam + triangleLoc1(p1) + p2;
225                                 int hoffset = at1 * numFreeParams + at2;
226
227                                 //mxLog("? H %d * g %d = p %d", hoffset, at2, at1);
228                                 matrixVectorProdTerm mvpt(hoffset, at2, at1);
229                                 state->hgProd.push_back(mvpt);
230
231                                 if (at1 != at2) {
232                                         matrixVectorProdTerm mvpt(hoffset, at1, at2);
233                                         state->hgProd.push_back(mvpt);
234                                 }
235
236                                 state->paramMap[at] = numFreeParams + hoffset;
237
238                                 state->ihessDivisor[at] =
239                                         state->paramLocations[at1] * state->paramLocations[at2];
240                         }
241                 }
242         }
243
244         state->haveItemMap = fc->varGroup->id;
245         //pia(state->paramMap.data(), state->itemDerivPadSize, itemParam->cols);
246 }
247
248 static double
249 ba81ComputeEMFit(omxFitFunction* oo, int want, FitContext *fc)
250 {
251         const double Scale = Global->llScale;
252         BA81FitState *state = (BA81FitState*) oo->argStruct;
253         BA81Expect *estate = (BA81Expect*) oo->expectation->argStruct;
254         omxMatrix *itemParam = estate->itemParam;
255         std::vector<const double*> &itemSpec = estate->itemSpec;
256         std::vector<int> &cumItemOutcomes = estate->cumItemOutcomes;
257         const int maxDims = estate->maxDims;
258         const size_t numItems = estate->itemSpec.size();
259         const int do_fit = want & FF_COMPUTE_FIT;
260         const int do_deriv = want & (FF_COMPUTE_GRADIENT | FF_COMPUTE_HESSIAN | FF_COMPUTE_IHESSIAN);
261
262         if (estate->verbose) mxLog("%s: em.fit(want fit=%d deriv=%d)", oo->matrix->name, do_fit, do_deriv);
263
264         if (do_fit) ba81OutcomeProb(estate, FALSE, TRUE);
265
266         const int thrDerivSize = itemParam->cols * state->itemDerivPadSize;
267         std::vector<double> thrDeriv(thrDerivSize * Global->numThreads);
268         double *wherePrep = estate->wherePrep.data();
269
270         double ll = 0;
271 #pragma omp parallel for num_threads(Global->numThreads) reduction(+:ll)
272         for (size_t ix=0; ix < numItems; ix++) {
273                 const int thrId = omx_absolute_thread_num();
274                 const double *spec = estate->itemSpec[ix];
275                 const int id = spec[RPF_ISpecID];
276                 const rpf_dLL1_t dLL1 = rpf_model[id].dLL1;
277                 const int iOutcomes = estate->itemOutcomes[ix];
278                 const int outcomeBase = cumItemOutcomes[ix] * estate->totalQuadPoints;
279                 const double *weight = estate->expected + outcomeBase;
280                 const double *oProb = estate->outcomeProb + outcomeBase;
281                 const double *iparam = omxMatrixColumn(itemParam, ix);
282                 double *myDeriv = thrDeriv.data() + thrDerivSize * thrId + ix * state->itemDerivPadSize;
283
284                 for (long qx=0; qx < estate->totalQuadPoints; qx++) {
285                         if (do_fit) {
286                                 for (int ox=0; ox < iOutcomes; ox++) {
287                                         ll += weight[ox] * oProb[ox];
288                                 }
289                         }
290                         if (do_deriv) {
291                                 (*dLL1)(spec, iparam, wherePrep + qx * maxDims, weight, myDeriv);
292                         }
293                         weight += iOutcomes;
294                         oProb += iOutcomes;
295                 }
296         }
297
298         size_t excluded = 0;
299
300         if (do_deriv) {
301                 double *deriv0 = thrDeriv.data();
302
303                 int perThread = itemParam->cols * state->itemDerivPadSize;
304                 for (int th=1; th < Global->numThreads; th++) {
305                         double *thrD = thrDeriv.data() + th * perThread;
306                         for (int ox=0; ox < perThread; ox++) deriv0[ox] += thrD[ox];
307                 }
308
309                 for (size_t ix=0; ix < numItems; ix++) {
310                         const double *spec = itemSpec[ix];
311                         int id = spec[RPF_ISpecID];
312                         double *iparam = omxMatrixColumn(itemParam, ix);
313                         double *pad = deriv0 + ix * state->itemDerivPadSize;
314                         (*rpf_model[id].dLL2)(spec, iparam, pad);
315                 }
316
317                 int numFreeParams = int(state->numFreeParam);
318                 int numParams = itemParam->cols * state->itemDerivPadSize;
319                 for (int ox=0; ox < numParams; ox++) {
320                         int to = state->paramMap[ox];
321                         if (to == -1) continue;
322
323                         // Need to check because this can happen if
324                         // lbounds/ubounds are not set appropriately.
325                         if (0 && !isfinite(deriv0[ox])) {
326                                 int item = ox / itemParam->rows;
327                                 mxLog("item parameters:\n");
328                                 const double *spec = itemSpec[item];
329                                 int id = spec[RPF_ISpecID];
330                                 int numParam = (*rpf_model[id].numParam)(spec);
331                                 double *iparam = omxMatrixColumn(itemParam, item);
332                                 pda(iparam, numParam, 1);
333                                 // Perhaps bounds can be pulled in from librpf? TODO
334                                 error("Deriv %d for item %d is %f; are you missing a lbound/ubound?",
335                                       ox, item, deriv0[ox]);
336                         }
337
338                         if (to < numFreeParams) {
339                                 if (want & FF_COMPUTE_GRADIENT) {
340                                         fc->grad[to] -= Scale * deriv0[ox];
341                                 }
342                         } else {
343                                 if (want & FF_COMPUTE_HESSIAN) {
344                                         int Hto = to - numFreeParams;
345                                         fc->hess[Hto] -= Scale * deriv0[ox];
346                                 }
347                         }
348                 }
349
350                 if (want & FF_COMPUTE_IHESSIAN) {
351                         for (size_t ix=0; ix < numItems; ix++) {
352                                 const double *spec = itemSpec[ix];
353                                 int id = spec[RPF_ISpecID];
354                                 int iParams = (*rpf_model[id].numParam)(spec);
355                                 double *pad = deriv0 + ix * state->itemDerivPadSize + iParams;
356                                 int *mask = state->itemParamFree.data() + ix * itemParam->rows;
357                                 double stress;
358                                 omxApproxInvertPackedPosDefTriangular(iParams, mask, pad, &stress);
359                                 // If items excluded then ihessDivisor is wrong TODO
360                                 if (stress) ++excluded;
361                         }
362                         for (int ox=0; ox < numParams; ox++) {
363                                 int to = state->paramMap[ox];
364                                 if (to == -1) continue;
365                                 if (to >= numFreeParams) {
366                                         int Hto = to - numFreeParams;
367                                         fc->ihess[Hto] -= deriv0[ox] / (Scale * state->ihessDivisor[ox]);
368                                 }
369                         }
370                 }
371         }
372
373         if (excluded && estate->verbose >= 1) {
374                 mxLog("%s: Hessian not positive definite for %lu/%lu items",
375                       oo->matrix->name, excluded, numItems);
376         }
377         if (excluded == numItems) {
378                 omxRaiseErrorf(globalState, "Hessian not positive definite for %lu/%lu items",
379                                excluded, numItems);
380         }
381
382         return Scale * ll;
383 }
384
385 void ba81SetFreeVarGroup(omxFitFunction *oo, FreeVarGroup *fvg)
386 {}
387
388 static bool approxInfo(omxFitFunction *oo, FitContext *fc)
389 {
390         const double abScale = fabs(Global->llScale);
391         omxExpectation *expectation = oo->expectation;
392         BA81FitState *state = (BA81FitState*) oo->argStruct;
393         BA81Expect *estate = (BA81Expect*) expectation->argStruct;
394         if (estate->verbose) mxLog("%s: approxInfo", oo->matrix->name);
395
396         ba81OutcomeProb(estate, FALSE, FALSE);
397
398         const int numThreads = Global->numThreads;
399         const int numUnique = estate->numUnique;
400         const int numSpecific = estate->numSpecific;
401         const int maxDims = estate->maxDims;
402         omxData *data = estate->data;
403         const int *rowMap = estate->rowMap;
404         int *numIdentical = estate->numIdentical;
405         const long totalQuadPoints = estate->totalQuadPoints;
406         omxMatrix *itemParam = estate->itemParam;
407         omxBuffer<double> patternLik(numUnique);
408
409         const int thrDerivSize = itemParam->cols * state->itemDerivPadSize;
410         const int totalOutcomes = estate->totalOutcomes;
411         const size_t numItems = estate->itemSpec.size();
412         const size_t numParam = fc->varGroup->vars.size();
413
414         if (numSpecific == 0) {
415                 omxBuffer<double> thrLxk(totalQuadPoints * numThreads);
416                 const double *wherePrep = estate->wherePrep.data();
417                 std::vector<double> thrBreadG(numThreads * numParam * numParam);
418                 std::vector<double> thrBreadH(numThreads * numParam * numParam);
419                 std::vector<double> thrMeat(numThreads * numParam * numParam);
420
421 #pragma omp parallel for num_threads(numThreads)
422                 for (int px=0; px < numUnique; px++) {
423                         int thrId = omx_absolute_thread_num();
424                         double *lxk = thrLxk.data() + thrId * totalQuadPoints;
425                         omxBuffer<double> deriv0(thrDerivSize);
426                         omxBuffer<double> expected(totalOutcomes); // can use maxOutcomes instead TODO
427                         double *breadG = thrBreadG.data() + thrId * numParam * numParam; //a
428                         double *breadH = thrBreadH.data() + thrId * numParam * numParam; //a
429                         double *meat = thrMeat.data() + thrId * numParam * numParam;   //b
430                         std::vector<double> patGrad(numParam);
431
432                         ba81LikelihoodSlow2(estate, px, lxk);
433
434                         // If patternLik is already valid, maybe could avoid this loop TODO
435                         double patternLik1 = 0;
436                         for (long qx=0; qx < totalQuadPoints; qx++) {
437                                 patternLik1 += lxk[qx];
438                         }
439                         patternLik[px] = patternLik1;
440
441                         // if (!validPatternLik(state, patternLik1))  complain
442
443                         double weight = 1 / patternLik[px];
444                         for (long qx=0; qx < totalQuadPoints; qx++) {
445                                 double tmp = lxk[qx] * weight;
446
447                                 OMXZERO(deriv0.data(), thrDerivSize);
448                                 std::vector<double> gradBuf(numParam);
449                                 int gradOffset = 0;
450
451                                 for (size_t ix=0; ix < numItems; ++ix) {
452                                         int pick = omxIntDataElementUnsafe(data, rowMap[px], ix);
453                                         if (pick == NA_INTEGER) continue;
454                                         pick -= 1;
455
456                                         const int iOutcomes = estate->itemOutcomes[ix];
457                                         OMXZERO(expected.data(), iOutcomes);
458                                         expected[pick] = 1;
459
460                                         const double *spec = estate->itemSpec[ix];
461                                         double *iparam = omxMatrixColumn(itemParam, ix);
462                                         const int id = spec[RPF_ISpecID];
463
464                                         double *myDeriv = deriv0.data() + ix * state->itemDerivPadSize;
465                                         (*rpf_model[id].dLL1)(spec, iparam, wherePrep + qx * maxDims,
466                                                               expected.data(), myDeriv);
467                                         (*rpf_model[id].dLL2)(spec, iparam, myDeriv);
468
469                                         for (int par = 0; par < state->paramPerItem[ix]; ++par) {
470                                                 int to = state->itemGradMap[gradOffset];
471                                                 if (to >= 0) gradBuf[to] -= myDeriv[par];
472                                                 ++gradOffset;
473                                         }
474                                 }
475
476                                 addSymOuterProd(abScale * tmp * numIdentical[px], gradBuf.data(), numParam, breadG);
477
478                                 for (size_t par=0; par < numParam; ++par) {
479                                         patGrad[par] += gradBuf[par] * tmp;
480                                 }
481                                 
482                                 for (int ox=0; ox < thrDerivSize; ox++) {
483                                         int to = state->paramMap[ox];
484                                         if (to == -1) continue;
485                                         if (to < int(numParam)) {
486                                                 // OK
487                                         } else {
488                                                 int Hto = to - numParam;
489                                                 breadH[Hto] += abScale * deriv0[ox] * tmp * numIdentical[px];
490                                         }
491                                 }
492                         }
493                         addSymOuterProd(abScale * numIdentical[px], patGrad.data(), numParam, meat);
494                 }
495
496                 // only need upper triangle TODO
497                 for (int tx=1; tx < numThreads; ++tx) {
498                         double *th = thrBreadG.data() + tx * numParam * numParam;
499                         for (size_t en=0; en < numParam * numParam; ++en) {
500                                 thrBreadG[en] += th[en];
501                         }
502                 }
503                 for (int tx=1; tx < numThreads; ++tx) {
504                         double *th = thrBreadH.data() + tx * numParam * numParam;
505                         for (size_t en=0; en < numParam * numParam; ++en) {
506                                 thrBreadH[en] += th[en];
507                         }
508                 }
509                 for (int tx=1; tx < numThreads; ++tx) {
510                         double *th = thrMeat.data() + tx * numParam * numParam;
511                         for (size_t en=0; en < numParam * numParam; ++en) {
512                                 thrMeat[en] += th[en];
513                         }
514                 }
515                 //pda(thrBreadG.data(), numParam, numParam);
516                 //pda(thrBreadH.data(), numParam, numParam);
517                 //pda(thrMeat.data(), numParam, numParam);
518                 if (fc->infoA) {
519                         for (size_t d1=0; d1 < numParam; ++d1) {
520                                 for (size_t d2=0; d2 < numParam; ++d2) {
521                                         int cell = d1 * numParam + d2;
522                                         fc->infoA[cell] -= thrMeat[cell] - (thrBreadG[cell] + thrBreadH[cell]);
523                                 }
524                         }
525                 }
526                 if (fc->infoB) {
527                         for (size_t d1=0; d1 < numParam; ++d1) {
528                                 for (size_t d2=0; d2 < numParam; ++d2) {
529                                         int cell = d1 * numParam + d2;
530                                         fc->infoB[cell] += thrMeat[cell];
531                                 }
532                         }
533                 }
534                 fc->sampleSize += data->rows; // remove? TODO
535         } else {
536                 const long totalPrimaryPoints = estate->totalPrimaryPoints;
537                 const long specificPoints = estate->quadGridSize;
538                 omxBuffer<double> thrLxk(totalQuadPoints * numSpecific * numThreads);
539                 omxBuffer<double> thrEi(totalPrimaryPoints * numThreads);
540                 omxBuffer<double> thrEis(totalPrimaryPoints * numSpecific * numThreads);
541
542 #pragma omp parallel for num_threads(numThreads)
543                 for (int px=0; px < numUnique; px++) {
544                         int thrId = omx_absolute_thread_num();
545                         double *lxk = thrLxk.data() + totalQuadPoints * numSpecific * thrId;
546                         double *Ei = thrEi.data() + totalPrimaryPoints * thrId;
547                         double *Eis = thrEis.data() + totalPrimaryPoints * numSpecific * thrId;
548                         cai2010EiEis(estate, px, lxk, Eis, Ei);
549
550                         for (long qloc=0, eisloc=0, qx=0; eisloc < totalPrimaryPoints * numSpecific; eisloc += numSpecific) {
551                                 for (long sx=0; sx < specificPoints; sx++) {
552                                         //double tmp0 = Eis[eisloc] * lxk[qloc];
553                                         // mapLatentDeriv(state, estate, tmp0,
554                                         //             derivCoef.data() + qx * derivPerPoint,
555                                         //             uniqueDeriv.data() + px * numLatents);
556
557                                         for (int Sgroup=0; Sgroup < numSpecific; Sgroup++) {
558                                                 //double lxk1 = lxk[qloc];
559                                                 //double Eis1 = Eis[eisloc + Sgroup];
560                                                 //double tmp = Eis1 * lxk1;
561                                                 // mapLatentDerivS(state, estate, Sgroup, tmp,
562                                                 //              derivCoef.data() + qx * derivPerPoint + priDerivCoef + 2 * Sgroup,
563                                                 //              uniqueDeriv.data() + px * numLatents);
564                                                 ++qloc;
565                                         }
566                                         ++qx;
567                                 }
568                         }
569
570                         double patternLik1 = 0;
571                         for (long qx=0; qx < totalPrimaryPoints; ++qx) {
572                                 patternLik1 += Ei[qx];
573                         }
574                         patternLik[px] = patternLik1;
575
576                         //double weight = numIdentical[px] / patternLik[px];
577                 }
578         }
579
580         return TRUE;
581 }
582
583 static void setLatentStartingValues(omxFitFunction *oo, FitContext *fc)
584 {
585         BA81FitState *state = (BA81FitState*) oo->argStruct;
586         BA81Expect *estate = (BA81Expect*) oo->expectation->argStruct;
587         std::vector<int> &latentMap = state->latentMap;
588         std::vector<double> &ElatentMean = estate->ElatentMean;
589         std::vector<double> &ElatentCov = estate->ElatentCov;
590         int maxAbilities = estate->maxAbilities;
591
592         if (!estate->Qpoint.size()) return; // if evaluating fit without estimating model
593         if (state->ElatentVersion == estate->ElatentVersion) return;
594
595         fc->changedEstimates = true;
596
597         for (int a1 = 0; a1 < maxAbilities; ++a1) {
598                 if (latentMap[a1] >= 0) {
599                         int to = latentMap[a1];
600                         fc->est[to] = ElatentMean[a1];
601                 }
602
603                 for (int a2 = 0; a2 <= a1; ++a2) {
604                         int to = latentMap[maxAbilities + triangleLoc1(a1) + a2];
605                         if (to < 0) continue;
606                         fc->est[to] = ElatentCov[a1 * maxAbilities + a2];
607                 }
608         }
609
610         state->ElatentVersion = estate->ElatentVersion;
611         //fc->log("setLatentStartingValues", FF_COMPUTE_ESTIMATE);
612 }
613
614 static void mapLatentDeriv(BA81FitState *state, BA81Expect *estate, double piece,
615                            double *derivCoef, double *derivOut)
616 {
617         const int maxAbilities = estate->maxAbilities;
618         const int pmax = estate->numSpecific? estate->maxDims - 1 : estate->maxDims;
619
620         int cx = 0;
621         for (int d1=0; d1 < pmax; ++d1) {
622                 double amt1 = piece * derivCoef[d1];
623                 derivOut[d1] += amt1;
624                 for (int d2=0; d2 <= d1; ++d2) {
625                         int to = maxAbilities + cx;
626                         double amt2 = piece * derivCoef[pmax + cx];
627                         derivOut[to] += amt2;
628                         ++cx;
629                 }
630         }
631 }
632
633 static void mapLatentDerivS(BA81FitState *state, BA81Expect *estate, int sgroup, double piece,
634                             double *derivCoef, double *derivOut)
635 {
636         int maxAbilities = estate->maxAbilities;
637         int maxDims = estate->maxDims;
638         int pmax = maxDims;
639         if (estate->numSpecific) pmax -= 1;
640
641         int sdim = pmax + sgroup;
642         double amt3 = piece * derivCoef[0];
643         derivOut[sdim] += amt3;
644
645         double amt4 = piece * derivCoef[1];
646         int to = maxAbilities + triangleLoc0(sdim);
647         derivOut[to] += amt4;
648 }
649
650 static void calcDerivCoef(BA81FitState *state, BA81Expect *estate, omxBuffer<double> &icov,
651                           const double *where, double *derivCoef)
652 {
653         omxMatrix *mean = estate->latentMeanOut;
654         omxMatrix *cov = estate->latentCovOut;
655         const int pDims = estate->numSpecific? estate->maxDims - 1 : estate->maxDims;
656         const char R='R';
657         const char L='L';
658         const char U='U';
659         const double alpha = 1;
660         const double beta = 0;
661         const int one = 1;
662
663         std::vector<double> whereDiff(pDims);
664         std::vector<double> whereGram(triangleLoc1(pDims));
665         for (int d1=0; d1 < pDims; ++d1) {
666                 whereDiff[d1] = where[d1] - omxVectorElement(mean, d1);
667         }
668         gramProduct(whereDiff.data(), whereDiff.size(), whereGram.data());
669
670         F77_CALL(dsymv)(&U, &pDims, &alpha, icov.data(), &pDims, whereDiff.data(), &one,
671                         &beta, derivCoef, &one);
672
673         std::vector<double> covGrad1(pDims * pDims);
674         std::vector<double> covGrad2(pDims * pDims);
675
676         int cx=0;
677         for (int d1=0; d1 < pDims; ++d1) {
678                 for (int d2=0; d2 <= d1; ++d2) {
679                         covGrad1[d2 * pDims + d1] = omxMatrixElement(cov, d2, d1) - whereGram[cx];
680                         ++cx;
681                 }
682         }
683
684         F77_CALL(dsymm)(&R, &L, &pDims, &pDims, &alpha, covGrad1.data(), &pDims, icov.data(),
685                         &pDims, &beta, covGrad2.data(), &pDims);
686         F77_CALL(dsymm)(&R, &L, &pDims, &pDims, &alpha, icov.data(), &pDims, covGrad2.data(),
687                         &pDims, &beta, covGrad1.data(), &pDims);
688
689         for (int d1=0; d1 < pDims; ++d1) {
690                 covGrad1[d1 * pDims + d1] /= 2.0;
691         }
692
693         cx = pDims;
694         for (int d1=0; d1 < pDims; ++d1) {
695                 int cell = d1 * pDims;
696                 for (int d2=0; d2 <= d1; ++d2) {
697                         derivCoef[cx] = -covGrad1[cell + d2];
698                         ++cx;
699                 }
700         }
701 }
702
703 static void calcDerivCoef1(BA81FitState *state, BA81Expect *estate,
704                            const double *where, int sgroup, double *derivCoef)
705 {
706         omxMatrix *mean = estate->latentMeanOut;
707         omxMatrix *cov = estate->latentCovOut;
708         const int maxDims = estate->maxDims;
709         const int specific = maxDims - 1 + sgroup;
710         double svar = omxMatrixElement(cov, specific, specific);
711         double whereDiff = where[maxDims-1] - omxVectorElement(mean, specific);
712         derivCoef[0] = whereDiff / svar;
713         derivCoef[1] = -(svar - whereDiff * whereDiff) / (2 * svar * svar);
714 }
715
716 static void xpd_finish_1pat(const double weight, const int numIdentical, const size_t numItems,
717                             const int numLatents, const size_t numParam,
718                             BA81FitState *state, BA81Expect *estate, omxMatrix *itemParam,
719                             std::vector<double> &deriv0, std::vector<double> &latentGrad,
720                             const double Scale,
721                             std::vector<double> &patGrad,
722                             double *grad, double *meat)
723 {
724         int gradOffset = 0;
725         for (size_t ix=0; ix < numItems; ++ix) {
726                 const double *spec = estate->itemSpec[ix];
727                 double *iparam = omxMatrixColumn(itemParam, ix);
728                 const int id = spec[RPF_ISpecID];
729                 double *myDeriv = deriv0.data() + ix * state->itemDerivPadSize;
730                 (*rpf_model[id].dLL2)(spec, iparam, myDeriv);
731
732                 for (int par = 0; par < state->paramPerItem[ix]; ++par) {
733                         int to = state->itemGradMap[gradOffset];
734                         if (to >= 0) patGrad[to] -= weight * myDeriv[par];
735                         ++gradOffset;
736                 }
737         }
738
739         for (int lx=0; lx < numLatents; ++lx) {
740                 int to = state->latentMap[lx];
741                 if (to >= 0) patGrad[to] += weight * latentGrad[lx];
742         }
743         for (size_t par=0; par < numParam; ++par) {
744                 grad[par] += patGrad[par] * Scale * numIdentical;
745         }
746         addSymOuterProd(fabs(Scale) * numIdentical, patGrad.data(), numParam, meat);
747 }
748
749 static bool xpd(omxFitFunction *oo, FitContext *fc)
750 {
751         const double Scale = Global->llScale;
752         omxExpectation *expectation = oo->expectation;
753         BA81FitState *state = (BA81FitState*) oo->argStruct;
754         BA81Expect *estate = (BA81Expect*) expectation->argStruct;
755         if (estate->verbose) mxLog("%s: cross product approximation", oo->matrix->name);
756
757         if (fc->infoMethod == INFO_METHOD_HESSIAN) {
758                 if (state->freeLatents) {
759                         omxRaiseErrorf(globalState, "Hessian not available with free latents");
760                         return FALSE;
761                 }
762
763                 ba81ComputeEMFit(oo, FF_COMPUTE_HESSIAN, fc);
764                 return TRUE;
765         }
766
767         if (fc->infoMethod != INFO_METHOD_MEAT) {
768                 omxRaiseErrorf(globalState, "Information matrix approximation method %d is not available",
769                                fc->infoMethod);
770                 return FALSE;
771         }
772
773         ba81OutcomeProb(estate, FALSE, FALSE);
774
775         const int numThreads = Global->numThreads;
776         const int numUnique = estate->numUnique;
777         const int numSpecific = estate->numSpecific;
778         const int maxDims = estate->maxDims;
779         const int pDims = numSpecific? maxDims-1 : maxDims;
780         const int maxAbilities = estate->maxAbilities;
781         omxMatrix *cov = estate->latentCovOut;
782         omxData *data = estate->data;
783         const int *rowMap = estate->rowMap;
784         int *numIdentical = estate->numIdentical;
785         const long totalQuadPoints = estate->totalQuadPoints;
786         omxMatrix *itemParam = estate->itemParam;
787         omxBuffer<double> patternLik(numUnique);
788
789         omxBuffer<double> icovBuffer(pDims * pDims);
790         for (int d1=0; d1 < pDims; ++d1) {
791                 for (int d2=0; d2 < pDims; ++d2) {
792                         icovBuffer[d1 * pDims + d2] = omxMatrixElement(cov, d1, d2);
793                 }
794         }
795         Matrix icovMat(icovBuffer.data(), pDims, pDims);
796         int info = InvertSymmetricPosDef(icovMat, 'U');
797         if (info) return FALSE;
798
799         // fill in rest from upper triangle
800         for (int rx=1; rx < pDims; ++rx) {
801                 for (int cx=0; cx < rx; ++cx) {
802                         icovBuffer[cx * pDims + rx] = icovBuffer[rx * pDims + cx];
803                 }
804         }
805
806         const int priDerivCoef = pDims + triangleLoc1(pDims);
807         const int numLatents = maxAbilities + triangleLoc1(maxAbilities);
808         const int thrDerivSize = itemParam->cols * state->itemDerivPadSize;
809         const int totalOutcomes = estate->totalOutcomes;
810         const size_t numItems = estate->itemSpec.size();
811         const size_t numParam = fc->varGroup->vars.size();
812         std::vector<double> thrGrad(numThreads * numParam);
813         std::vector<double> thrMeat(numThreads * numParam * numParam);
814         const double *wherePrep = estate->wherePrep.data();
815
816         if (numSpecific == 0) {
817                 omxBuffer<double> thrLxk(totalQuadPoints * numThreads);
818                 omxBuffer<double> derivCoef(totalQuadPoints * priDerivCoef);
819
820 #pragma omp parallel for num_threads(numThreads)
821                 for (long qx=0; qx < totalQuadPoints; qx++) {
822                         const double *where = wherePrep + qx * maxDims;
823                         calcDerivCoef(state, estate, icovBuffer, where,
824                                       derivCoef.data() + qx * priDerivCoef);
825                 }
826
827 #pragma omp parallel for num_threads(numThreads)
828                 for (int px=0; px < numUnique; px++) {
829                         int thrId = omx_absolute_thread_num();
830                         double *lxk = thrLxk.data() + thrId * totalQuadPoints;
831                         omxBuffer<double> expected(totalOutcomes); // can use maxOutcomes instead TODO
832                         std::vector<double> deriv0(thrDerivSize);
833                         std::vector<double> latentGrad(numLatents);
834                         std::vector<double> patGrad(numParam);
835                         double *grad = thrGrad.data() + thrId * numParam;
836                         double *meat = thrMeat.data() + thrId * numParam * numParam;
837                         ba81LikelihoodSlow2(estate, px, lxk);
838
839                         // If patternLik is already valid, maybe could avoid this loop TODO
840                         double patternLik1 = 0;
841                         for (long qx=0; qx < totalQuadPoints; qx++) {
842                                 patternLik1 += lxk[qx];
843                         }
844                         patternLik[px] = patternLik1;
845
846                         // if (!validPatternLik(state, patternLik1))  complain, TODO
847
848                         for (long qx=0; qx < totalQuadPoints; qx++) {
849                                 double tmp = lxk[qx];
850                                 mapLatentDeriv(state, estate, tmp, derivCoef.data() + qx * priDerivCoef,
851                                                latentGrad.data());
852
853                                 for (size_t ix=0; ix < numItems; ++ix) {
854                                         int pick = omxIntDataElementUnsafe(data, rowMap[px], ix);
855                                         if (pick == NA_INTEGER) continue;
856                                         OMXZERO(expected.data(), estate->itemOutcomes[ix]);
857                                         expected[pick-1] = tmp;
858                                         const double *spec = estate->itemSpec[ix];
859                                         double *iparam = omxMatrixColumn(itemParam, ix);
860                                         const int id = spec[RPF_ISpecID];
861                                         double *myDeriv = deriv0.data() + ix * state->itemDerivPadSize;
862                                         (*rpf_model[id].dLL1)(spec, iparam, wherePrep + qx * maxDims,
863                                                               expected.data(), myDeriv);
864                                 }
865                         }
866
867                         xpd_finish_1pat(1 / patternLik1, numIdentical[px], numItems, numLatents, numParam,
868                                         state, estate, itemParam, deriv0, latentGrad, Scale, patGrad, grad, meat);
869                 }
870         } else {
871                 const long totalPrimaryPoints = estate->totalPrimaryPoints;
872                 const long specificPoints = estate->quadGridSize;
873                 omxBuffer<double> thrLxk(totalQuadPoints * numSpecific * numThreads);
874                 omxBuffer<double> thrEi(totalPrimaryPoints * numThreads);
875                 omxBuffer<double> thrEis(totalPrimaryPoints * numSpecific * numThreads);
876                 const int derivPerPoint = priDerivCoef + 2 * numSpecific;
877                 omxBuffer<double> derivCoef(totalQuadPoints * derivPerPoint);
878
879 #pragma omp parallel for num_threads(numThreads)
880                 for (long qx=0; qx < totalQuadPoints; qx++) {
881                         const double *where = wherePrep + qx * maxDims;
882                         calcDerivCoef(state, estate, icovBuffer, where,
883                                       derivCoef.data() + qx * derivPerPoint);
884                         for (int Sgroup=0; Sgroup < numSpecific; ++Sgroup) {
885                                 calcDerivCoef1(state, estate, where, Sgroup,
886                                                derivCoef.data() + qx * derivPerPoint + priDerivCoef + 2 * Sgroup);
887                         }
888                 }
889
890 #pragma omp parallel for num_threads(numThreads)
891                 for (int px=0; px < numUnique; px++) {
892                         int thrId = omx_absolute_thread_num();
893                         double *lxk = thrLxk.data() + totalQuadPoints * numSpecific * thrId;
894                         double *Ei = thrEi.data() + totalPrimaryPoints * thrId;
895                         double *Eis = thrEis.data() + totalPrimaryPoints * numSpecific * thrId;
896                         omxBuffer<double> expected(totalOutcomes); // can use maxOutcomes instead TODO
897                         std::vector<double> deriv0(thrDerivSize);
898                         std::vector<double> latentGrad(numLatents);
899                         std::vector<double> patGrad(numParam);
900                         double *grad = thrGrad.data() + thrId * numParam;
901                         double *meat = thrMeat.data() + thrId * numParam * numParam;
902                         cai2010EiEis(estate, px, lxk, Eis, Ei);
903
904                         for (long qloc=0, eisloc=0, qx=0; eisloc < totalPrimaryPoints * numSpecific; eisloc += numSpecific) {
905                                 for (long sx=0; sx < specificPoints; sx++) {
906                                         mapLatentDeriv(state, estate, Eis[eisloc] * lxk[qloc],
907                                                        derivCoef.data() + qx * derivPerPoint,
908                                                        latentGrad.data());
909
910                                         for (int Sgroup=0; Sgroup < numSpecific; Sgroup++) {
911                                                 double lxk1 = lxk[qloc];
912                                                 double Eis1 = Eis[eisloc + Sgroup];
913                                                 double tmp = Eis1 * lxk1;
914                                                 mapLatentDerivS(state, estate, Sgroup, tmp,
915                                                                 derivCoef.data() + qx * derivPerPoint + priDerivCoef + 2 * Sgroup,
916                                                                 latentGrad.data());
917
918                                                 for (size_t ix=0; ix < numItems; ++ix) {
919                                                         if (estate->Sgroup[ix] != Sgroup) continue;
920                                                         int pick = omxIntDataElementUnsafe(data, rowMap[px], ix);
921                                                         if (pick == NA_INTEGER) continue;
922                                                         OMXZERO(expected.data(), estate->itemOutcomes[ix]);
923                                                         expected[pick-1] = tmp;
924                                                         const double *spec = estate->itemSpec[ix];
925                                                         double *iparam = omxMatrixColumn(itemParam, ix);
926                                                         const int id = spec[RPF_ISpecID];
927                                                         double *myDeriv = deriv0.data() + ix * state->itemDerivPadSize;
928                                                         (*rpf_model[id].dLL1)(spec, iparam, wherePrep + qx * maxDims,
929                                                                               expected.data(), myDeriv);
930                                                 }
931                                                 ++qloc;
932                                         }
933                                         ++qx;
934                                 }
935                         }
936
937                         // If patternLik is already valid, maybe could avoid this loop TODO
938                         double patternLik1 = 0;
939                         for (long qx=0; qx < totalPrimaryPoints; ++qx) {
940                                 patternLik1 += Ei[qx];
941                         }
942                         patternLik[px] = patternLik1;
943
944                         xpd_finish_1pat(1 / patternLik1, numIdentical[px], numItems, numLatents, numParam,
945                                         state, estate, itemParam, deriv0, latentGrad, Scale, patGrad, grad, meat);
946                 }
947         }
948
949         for (int tx=1; tx < numThreads; ++tx) {
950                 double *th = thrGrad.data() + tx * numParam;
951                 for (size_t en=0; en < numParam; ++en) {
952                         thrGrad[en] += th[en];
953                 }
954         }
955         for (int tx=1; tx < numThreads; ++tx) {
956                 double *th = thrMeat.data() + tx * numParam * numParam;
957                 for (size_t en=0; en < numParam * numParam; ++en) {
958                         thrMeat[en] += th[en];
959                 }
960         }
961         for (size_t d1=0; d1 < numParam; ++d1) {
962                 fc->grad[d1] += thrGrad[d1];
963                 for (size_t d2=0; d2 < numParam; ++d2) {
964                         int cell = d1 * numParam + d2;
965                         fc->infoB[cell] += thrMeat[cell];
966                 }
967         }
968
969         return TRUE;
970 }
971
972 static double
973 ba81ComputeFit(omxFitFunction* oo, int want, FitContext *fc)
974 {
975         BA81FitState *state = (BA81FitState*) oo->argStruct;
976         BA81Expect *estate = (BA81Expect*) oo->expectation->argStruct;
977
978         if (estate->type == EXPECTATION_AUGMENTED) {
979                 buildItemParamMap(oo, fc);
980
981                 if (want & FF_COMPUTE_PARAMFLAVOR) {
982                         for (size_t px=0; px < state->numFreeParam; ++px) {
983                                 if (state->paramFlavor[px] < 0) continue;
984                                 fc->flavor[px] = state->paramFlavor[px];
985                         }
986                         return 0;
987                 }
988
989                 if (want & FF_COMPUTE_HGPROD) {
990                         for (size_t px=0; px < state->hgProd.size(); ++px) {
991                                 fc->hgProd.push_back(state->hgProd[px]);
992                         }
993                         return 0;
994                 }
995
996                 if (want & FF_COMPUTE_PREOPTIMIZE) {
997                         omxExpectationCompute(oo->expectation, NULL);
998                         // schilling_bock_2005_rescale(oo, fc); seems counterproductive
999                         return 0;
1000                 }
1001
1002                 if (want & FF_COMPUTE_INFO) {
1003                         buildLatentParamMap(oo, fc);
1004                         buildItemParamMap(oo, fc);
1005                         ba81SetupQuadrature(oo->expectation);
1006                         if (!xpd(oo, fc)) {
1007                                 return INFINITY;
1008                         }
1009                         return 0;
1010                 }
1011
1012                 double got = ba81ComputeEMFit(oo, want, fc);
1013                 return got;
1014         } else if (estate->type == EXPECTATION_OBSERVED) {
1015
1016                 if (want & FF_COMPUTE_PREOPTIMIZE) {
1017                         buildLatentParamMap(oo, fc);
1018                         if (state->freeLatents) {
1019                                 setLatentStartingValues(oo, fc);
1020                         }
1021                         return 0;
1022                 }
1023
1024                 if (want & FF_COMPUTE_INFO) {
1025                         buildLatentParamMap(oo, fc); // only to check state->freeLatents
1026                         buildItemParamMap(oo, fc);
1027
1028                         if (state->freeLatents) {
1029                                 omxRaiseErrorf(globalState, "Cannot approximate latent parameter gradients");
1030                         } else {
1031                                 ba81SetupQuadrature(oo->expectation);
1032                                 if (!approxInfo(oo, fc)) {
1033                                         return INFINITY;
1034                                 }
1035                         }
1036                 }
1037                 if (want & FF_COMPUTE_HESSIAN) {
1038                         warning("%s: Hessian not available", oo->matrix->name);
1039                 }
1040
1041                 if (want & FF_COMPUTE_MAXABSCHANGE) {
1042                         double mac = std::max(omxMaxAbsDiff(state->itemParam, estate->itemParam),
1043                                               omxMaxAbsDiff(state->latentMean, estate->latentMeanOut));
1044                         fc->mac = std::max(mac, omxMaxAbsDiff(state->latentCov, estate->latentCovOut));
1045                         state->copyEstimates(estate);
1046                 }
1047
1048                 if (want & FF_COMPUTE_FIT) {
1049                         omxExpectationCompute(oo->expectation, NULL);
1050
1051                         double *patternLik = estate->patternLik;
1052                         int *numIdentical = estate->numIdentical;
1053                         int numUnique = estate->numUnique;
1054                         estate->excludedPatterns = 0;
1055                         const double LogLargest = estate->LogLargestDouble;
1056                         double got = 0;
1057 #pragma omp parallel for num_threads(Global->numThreads) reduction(+:got)
1058                         for (int ux=0; ux < numUnique; ux++) {
1059                                 if (!validPatternLik(estate, patternLik[ux])) {
1060 #pragma omp atomic
1061                                         ++estate->excludedPatterns;
1062                                         // somehow indicate that this -2LL is provisional TODO
1063                                         continue;
1064                                 }
1065                                 got += numIdentical[ux] * (log(patternLik[ux]) - LogLargest);
1066                         }
1067                         if (estate->verbose) mxLog("%s: fit (%d/%d excluded)",
1068                                                    oo->matrix->name, estate->excludedPatterns, numUnique);
1069                         //mxLog("fit %.4f", -2 * got);
1070                         return Global->llScale * got;
1071                 }
1072
1073                 return 0;
1074         } else {
1075                 error("%s: Select EM or regular mode before computing %d", oo->matrix->name, want);
1076         }
1077 }
1078
1079 static void ba81Compute(omxFitFunction *oo, int want, FitContext *fc)
1080 {
1081         if (!want) return;
1082         double got = ba81ComputeFit(oo, want, fc);
1083         if (got) oo->matrix->data[0] = got;
1084 }
1085
1086 BA81FitState::~BA81FitState()
1087 {
1088         omxFreeAllMatrixData(itemParam);
1089         omxFreeAllMatrixData(latentMean);
1090         omxFreeAllMatrixData(latentCov);
1091 }
1092
1093 static void ba81Destroy(omxFitFunction *oo) {
1094         BA81FitState *state = (BA81FitState *) oo->argStruct;
1095         delete state;
1096 }
1097
1098 void omxInitFitFunctionBA81(omxFitFunction* oo)
1099 {
1100         if (!oo->argStruct) { // ugh!
1101                 BA81FitState *state = new BA81FitState;
1102                 oo->argStruct = state;
1103         }
1104
1105         BA81FitState *state = (BA81FitState*) oo->argStruct;
1106
1107         omxExpectation *expectation = oo->expectation;
1108         BA81Expect *estate = (BA81Expect*) expectation->argStruct;
1109
1110         //newObj->data = oo->expectation->data;
1111
1112         oo->computeFun = ba81Compute;
1113         oo->setVarGroup = ba81SetFreeVarGroup;
1114         oo->destructFun = ba81Destroy;
1115         oo->gradientAvailable = TRUE;
1116         oo->hessianAvailable = TRUE;
1117         oo->parametersHaveFlavor = TRUE;
1118
1119         int maxParam = estate->itemParam->rows;
1120         state->itemDerivPadSize = maxParam + triangleLoc1(maxParam);
1121
1122         int numItems = estate->itemParam->cols;
1123         for (int ix=0; ix < numItems; ix++) {
1124                 const double *spec = estate->itemSpec[ix];
1125                 int id = spec[RPF_ISpecID];
1126                 if (id < 0 || id >= rpf_numModels) {
1127                         error("ItemSpec %d has unknown item model %d", ix, id);
1128                 }
1129         }
1130
1131         state->itemParam = omxInitMatrix(NULL, 0, 0, TRUE, globalState);
1132         state->latentMean = omxInitMatrix(NULL, 0, 0, TRUE, globalState);
1133         state->latentCov = omxInitMatrix(NULL, 0, 0, TRUE, globalState);
1134         state->copyEstimates(estate);
1135         state->ElatentVersion = -1;
1136 }