2 Copyright 2012-2013 Joshua Nathaniel Pritikin and contributors
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.
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.
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/>.
18 #include "omxFitFunction.h"
19 #include "omxExpectationBA81.h"
20 #include "omxOpenmpWrap.h"
21 #include "libifa-rpf.h"
23 static const char *NAME = "FitFunctionBA81";
27 omxMatrix *itemParam; // M step version
28 int derivPadSize; // maxParam + maxParam*(1+maxParam)/2
29 double *thrDeriv; // itemParam->cols * derivPadSize * thread
30 int *paramMap; // itemParam->cols * derivPadSize -> index of free parameter
31 std::vector<int> latentMeanMap;
32 std::vector<int> latentCovMap;
33 std::vector<int> NAtriangle;
35 omxMatrix *customPrior;
37 double *tmpLatentMean; // maxDims
38 double *tmpLatentCov; // maxDims * maxDims ; only lower triangle is used
42 std::vector< FreeVarGroup* > varGroups;
43 FreeVarGroup *latentFVG;
48 BA81FitState::BA81FitState()
61 static void buildParamMap(omxFitFunction* oo)
63 BA81FitState *state = (BA81FitState *) oo->argStruct;
64 BA81Expect *estate = (BA81Expect*) oo->expectation->argStruct;
65 omxMatrix *itemParam = state->itemParam;
66 int size = itemParam->cols * state->derivPadSize;
67 int maxAbilities = estate->maxAbilities;
68 int meanNum = estate->latentMeanOut->matrixNumber;
69 int covNum = estate->latentCovOut->matrixNumber;
70 FreeVarGroup *fvg = state->latentFVG;
72 state->latentMeanMap.assign(maxAbilities, -1);
73 state->latentCovMap.assign(maxAbilities * maxAbilities, -1);
75 for (size_t px=0; px < fvg->vars.size(); px++) {
76 omxFreeVar *fv = fvg->vars[px];
77 for (size_t lx=0; lx < fv->locations.size(); lx++) {
78 omxFreeVarLocation *loc = &fv->locations[lx];
79 int matNum = ~loc->matrix;
80 if (matNum == meanNum) {
81 state->latentMeanMap[loc->row + loc->col] = px;
82 } else if (matNum == covNum) {
83 state->latentCovMap[loc->col * maxAbilities + loc->row] = px;
88 state->paramMap = Realloc(NULL, size, int); // matrix location to free param index
89 for (int px=0; px < size; px++) {
90 state->paramMap[px] = -1;
93 size_t numFreeParams = oo->freeVarGroup->vars.size();
94 int *pRow = Realloc(NULL, numFreeParams, int);
95 int *pCol = Realloc(NULL, numFreeParams, int);
97 for (size_t px=0; px < numFreeParams; px++) {
100 omxFreeVar *fv = oo->freeVarGroup->vars[px];
101 for (size_t lx=0; lx < fv->locations.size(); lx++) {
102 omxFreeVarLocation *loc = &fv->locations[lx];
103 int matNum = ~loc->matrix;
104 if (matNum == itemParam->matrixNumber) {
107 int at = pCol[px] * state->derivPadSize + pRow[px];
108 state->paramMap[at] = px;
110 const double *spec = omxMatrixColumn(estate->itemSpec, loc->col);
111 int id = spec[RPF_ISpecID];
113 (*rpf_model[id].paramBound)(spec, loc->row, &upper, &lower);
114 if (fv->lbound == NEG_INF && isfinite(lower)) fv->lbound = lower;
115 if (fv->ubound == INF && isfinite(upper)) fv->ubound = upper;
120 for (size_t p1=0; p1 < numFreeParams; p1++) {
121 for (size_t p2=p1; p2 < numFreeParams; p2++) {
122 if (pCol[p1] == -1 || pCol[p1] != pCol[p2]) continue;
123 const double *spec = omxMatrixColumn(estate->itemSpec, pCol[p1]);
124 int id = spec[RPF_ISpecID];
125 int numParam = (*rpf_model[id].numParam)(spec);
128 if (r1 > r2) { int tmp=r1; r1=r2; r2=tmp; }
130 for (int rx=1; rx <= r2; rx++) rowOffset += rx;
131 int at = pCol[p1] * state->derivPadSize + numParam + rowOffset + r1;
132 state->paramMap[at] = numFreeParams + p1 * numFreeParams + p2;
133 if (p2 != p1) state->NAtriangle.push_back(p2 * numFreeParams + p1);
140 state->thrDeriv = Realloc(NULL, itemParam->cols * state->derivPadSize * Global->numThreads, double);
143 OMXINLINE static double
144 ba81Fit1Ordinate(omxFitFunction* oo, const int *quad, const double *weight, int want)
146 BA81FitState *state = (BA81FitState*) oo->argStruct;
147 BA81Expect *estate = (BA81Expect*) oo->expectation->argStruct;
148 omxMatrix *itemSpec = estate->itemSpec;
149 omxMatrix *itemParam = state->itemParam;
150 int numItems = itemParam->cols;
151 int maxOutcomes = estate->maxOutcomes;
152 int maxDims = estate->maxDims;
153 double *myDeriv = state->thrDeriv + itemParam->cols * state->derivPadSize * omx_absolute_thread_num();
154 int do_deriv = want & (FF_COMPUTE_GRADIENT | FF_COMPUTE_HESSIAN);
156 double where[maxDims];
157 pointToWhere(estate, quad, where, maxDims);
159 double *outcomeProb = computeRPF(estate, itemParam, quad); // avoid malloc/free? TODO
160 if (!outcomeProb) return 0;
163 for (int ix=0; ix < numItems; ix++) {
164 const double *spec = omxMatrixColumn(itemSpec, ix);
165 int id = spec[RPF_ISpecID];
166 int iOutcomes = spec[RPF_ISpecOutcomes];
168 double area = exp(logAreaProduct(estate, quad, estate->Sgroup[ix])); // avoid exp() here? TODO
169 for (int ox=0; ox < iOutcomes; ox++) {
171 #pragma omp critical(ba81Fit1OrdinateDebug1)
172 if (!std::isfinite(outcomeProb[ix * maxOutcomes + ox])) {
173 pda(itemParam->data, itemParam->rows, itemParam->cols);
174 pda(outcomeProb, outcomes, numItems);
175 error("RPF produced NAs");
178 double got = weight[ox] * outcomeProb[ix * maxOutcomes + ox];
179 thr_ll += got * area;
183 double *iparam = omxMatrixColumn(itemParam, ix);
184 double *pad = myDeriv + ix * state->derivPadSize;
185 (*rpf_model[id].dLL1)(spec, iparam, where, area, weight, pad);
196 ba81ComputeMFit1(omxFitFunction* oo, int want, double *gradient, double *hessian)
198 BA81FitState *state = (BA81FitState*) oo->argStruct;
199 BA81Expect *estate = (BA81Expect*) oo->expectation->argStruct;
200 omxMatrix *customPrior = state->customPrior;
201 omxMatrix *itemParam = state->itemParam;
202 omxMatrix *itemSpec = estate->itemSpec;
203 int maxDims = estate->maxDims;
204 const int totalOutcomes = estate->totalOutcomes;
208 omxRecompute(customPrior);
209 ll = customPrior->data[0];
210 // need deriv adjustment TODO
214 omxPrint(itemParam, "item param");
215 error("Bayesian prior returned %g; do you need to add a lbound/ubound?", ll);
218 #pragma omp parallel for num_threads(Global->numThreads)
219 for (long qx=0; qx < estate->totalQuadPoints; qx++) {
220 //double area = exp(state->priLogQarea[qx]); // avoid exp() here? TODO
222 decodeLocation(qx, maxDims, estate->quadGridSize, quad);
223 double *weight = estate->expected + qx * totalOutcomes;
224 double thr_ll = ba81Fit1Ordinate(oo, quad, weight, want);
231 double *deriv0 = state->thrDeriv;
233 int perThread = itemParam->cols * state->derivPadSize;
234 for (int th=1; th < Global->numThreads; th++) {
235 double *thrD = state->thrDeriv + th * perThread;
236 for (int ox=0; ox < perThread; ox++) deriv0[ox] += thrD[ox];
239 int numItems = itemParam->cols;
240 for (int ix=0; ix < numItems; ix++) {
241 const double *spec = omxMatrixColumn(itemSpec, ix);
242 int id = spec[RPF_ISpecID];
243 double *iparam = omxMatrixColumn(itemParam, ix);
244 double *pad = deriv0 + ix * state->derivPadSize;
245 (*rpf_model[id].dLL2)(spec, iparam, pad);
248 int numFreeParams = int(oo->freeVarGroup->vars.size());
249 int numParams = itemParam->cols * state->derivPadSize;
250 for (int ox=0; ox < numParams; ox++) {
251 int to = state->paramMap[ox];
252 if (to == -1) continue;
254 // Need to check because this can happen if
255 // lbounds/ubounds are not set appropriately.
256 if (0 && !isfinite(deriv0[ox])) {
257 int item = ox / itemParam->rows;
258 mxLog("item parameters:\n");
259 const double *spec = omxMatrixColumn(itemSpec, item);
260 int id = spec[RPF_ISpecID];
261 int numParam = (*rpf_model[id].numParam)(spec);
262 double *iparam = omxMatrixColumn(itemParam, item);
263 pda(iparam, numParam, 1);
264 // Perhaps bounds can be pulled in from librpf? TODO
265 error("Deriv %d for item %d is %f; are you missing a lbound/ubound?",
266 ox, item, deriv0[ox]);
269 if (to < numFreeParams) {
270 gradient[to] -= deriv0[ox];
272 hessian[to - numFreeParams] -= deriv0[ox];
281 moveLatentDistribution(omxFitFunction *oo, FitContext *fc,
282 double *ElatentMean, double *ElatentCov)
284 BA81FitState *state = (BA81FitState*) oo->argStruct;
285 BA81Expect *estate = (BA81Expect*) oo->expectation->argStruct;
286 omxMatrix *itemSpec = estate->itemSpec;
287 omxMatrix *itemParam = state->itemParam;
288 omxMatrix *design = estate->design;
289 double *tmpLatentMean = state->tmpLatentMean;
290 double *tmpLatentCov = state->tmpLatentCov;
291 int maxDims = estate->maxDims;
292 int maxAbilities = estate->maxAbilities;
294 int numItems = itemParam->cols;
295 for (int ix=0; ix < numItems; ix++) {
296 const double *spec = omxMatrixColumn(itemSpec, ix);
297 int id = spec[RPF_ISpecID];
298 const double *rawDesign = omxMatrixColumn(design, ix);
299 int idesign[design->rows];
301 for (int dx=0; dx < design->rows; dx++) {
302 if (isfinite(rawDesign[dx])) {
303 idesign[idx++] = rawDesign[dx]-1;
308 for (int d1=0; d1 < idx; d1++) {
309 if (idesign[d1] == -1) {
310 tmpLatentMean[d1] = 0;
312 tmpLatentMean[d1] = ElatentMean[idesign[d1]];
314 for (int d2=0; d2 <= d1; d2++) {
315 int cell = idesign[d2] * maxAbilities + idesign[d1];
316 if (idesign[d1] == -1 || idesign[d2] == -1) {
317 tmpLatentCov[d2 * maxDims + d1] = d1==d2? 1 : 0;
319 tmpLatentCov[d2 * maxDims + d1] = ElatentCov[cell];
323 if (1) { // ease debugging, make optional TODO
324 for (int d1=idx; d1 < maxDims; d1++) tmpLatentMean[d1] = nan("");
325 for (int d1=0; d1 < maxDims; d1++) {
326 for (int d2=0; d2 < maxDims; d2++) {
327 if (d1 < idx && d2 < idx) continue;
328 tmpLatentCov[d2 * maxDims + d1] = nan("");
332 double *iparam = omxMatrixColumn(itemParam, ix);
333 int *mask = state->paramMap + state->derivPadSize * ix;
334 rpf_model[id].rescale(spec, iparam, mask, tmpLatentMean, tmpLatentCov);
337 int numFreeParams = int(oo->freeVarGroup->vars.size());
338 for (int rx=0; rx < itemParam->rows; rx++) {
339 for (int cx=0; cx < itemParam->cols; cx++) {
340 int vx = state->paramMap[cx * state->derivPadSize + rx];
341 if (vx >= 0 && vx < numFreeParams) {
342 fc->est[vx] = omxMatrixElement(itemParam, rx, cx);
349 schilling_bock_2005_rescale(omxFitFunction *oo, FitContext *fc)
351 BA81FitState *state = (BA81FitState*) oo->argStruct;
352 BA81Expect *estate = (BA81Expect*) oo->expectation->argStruct;
353 double *ElatentMean = estate->ElatentMean;
354 double *ElatentCov = estate->ElatentCov;
355 int maxAbilities = estate->maxAbilities;
357 //mxLog("schilling bock\n");
358 //pda(ElatentMean, maxAbilities, 1);
359 //pda(ElatentCov, maxAbilities, maxAbilities);
360 //omxPrint(design, "design");
362 // use omxDPOTRF instead? TODO
363 const char triangle = 'L';
364 F77_CALL(dpotrf)(&triangle, &maxAbilities, ElatentCov, &maxAbilities, &state->choleskyError);
365 if (state->choleskyError != 0) {
366 warning("Cholesky failed with %d; rescaling disabled", state->choleskyError); // make error TODO?
370 moveLatentDistribution(oo, fc, ElatentMean, ElatentCov);
371 fc->copyParamToModel(globalState);
374 OMXINLINE static void
375 updateLatentParam(omxFitFunction* oo, FitContext *fc)
377 BA81FitState *state = (BA81FitState*) oo->argStruct;
378 BA81Expect *estate = (BA81Expect*) oo->expectation->argStruct;
379 int maxAbilities = estate->maxAbilities;
381 // TODO need denom for multigroup
382 for (int a1=0; a1 < maxAbilities; ++a1) {
383 if (state->latentMeanMap[a1] >= 0) {
384 double val = estate->ElatentMean[a1];
385 int vx = state->latentMeanMap[a1];
386 omxFreeVar *fv = fc->varGroup->vars[vx];
387 if (val < fv->lbound) val = fv->lbound;
388 if (val > fv->ubound) val = fv->ubound;
391 for (int a2=0; a2 < maxAbilities; ++a2) {
392 int cell = a2 * maxAbilities + a1;
393 if (state->latentCovMap[cell] < 0) continue;
394 double val = estate->ElatentCov[cell];
395 int vx = state->latentCovMap[cell];
396 omxFreeVar *fv = fc->varGroup->vars[vx];
397 if (val < fv->lbound) val = fv->lbound;
398 if (val > fv->ubound) val = fv->ubound;
403 fc->copyParamToModel(globalState);
406 void ba81SetFreeVarGroup(omxFitFunction *oo, FreeVarGroup *fvg) // too ad hoc? TODO
408 if (!oo->argStruct) { // ugh!
409 BA81FitState *state = new BA81FitState;
410 oo->argStruct = state;
413 BA81FitState *state = (BA81FitState*) oo->argStruct;
415 state->varGroups.push_back(fvg);
416 if (state->varGroups.size() == 2) {
418 if (state->varGroups[0] == state->varGroups[1])
419 warning("Cannot recognize correct free parameter groups");
420 if (state->varGroups[0]->vars.size() > state->varGroups[1]->vars.size())
422 oo->freeVarGroup = state->varGroups[small];
423 state->latentFVG = state->varGroups[!small];
424 } else if (state->varGroups.size() > 2) {
430 ba81ComputeFit(omxFitFunction* oo, int want, FitContext *fc)
432 BA81FitState *state = (BA81FitState*) oo->argStruct;
434 if (!state->paramMap) buildParamMap(oo);
436 if (want & FF_COMPUTE_PREOPTIMIZE) {
437 if (state->rescale) schilling_bock_2005_rescale(oo, fc); // how does this work in multigroup? TODO
441 BA81Expect *estate = (BA81Expect*) oo->expectation->argStruct;
445 if (estate->type == EXPECTATION_AUGMENTED) {
446 if (fc->varGroup != oo->freeVarGroup) error("FreeVarGroup mismatch");
448 if (want & FF_COMPUTE_GRADIENT) ++state->gradientCount;
450 omxMatrix *itemParam = state->itemParam;
451 OMXZERO(state->thrDeriv, state->derivPadSize * itemParam->cols * Global->numThreads);
453 for (size_t nx=0; nx < state->NAtriangle.size(); ++nx) {
454 fc->hess[ state->NAtriangle[nx] ] = nan("symmetric");
457 double got = ba81ComputeMFit1(oo, want, fc->grad, fc->hess);
459 } else if (estate->type == EXPECTATION_OBSERVED) {
460 if (fc->varGroup != state->latentFVG) error("FreeVarGroup mismatch");
462 updateLatentParam(oo, fc);
464 double *patternLik = estate->patternLik;
465 int *numIdentical = estate->numIdentical;
466 int numUnique = estate->numUnique;
468 for (int ux=0; ux < numUnique; ux++) {
469 got += numIdentical[ux] * patternLik[ux];
475 static void ba81Compute(omxFitFunction *oo, int want, FitContext *fc)
478 oo->matrix->data[0] = ba81ComputeFit(oo, want, fc);
481 static void ba81Destroy(omxFitFunction *oo) {
482 BA81FitState *state = (BA81FitState *) oo->argStruct;
484 omxFreeAllMatrixData(state->customPrior);
485 Free(state->paramMap);
486 Free(state->thrDeriv);
487 Free(state->tmpLatentMean);
488 Free(state->tmpLatentCov);
489 omxFreeAllMatrixData(state->itemParam);
493 void omxInitFitFunctionBA81(omxFitFunction* oo)
495 if (!oo->argStruct) { // ugh!
496 BA81FitState *state = new BA81FitState;
497 oo->argStruct = state;
500 BA81FitState *state = (BA81FitState*) oo->argStruct;
501 SEXP rObj = oo->rObj;
503 omxExpectation *expectation = oo->expectation;
504 BA81Expect *estate = (BA81Expect*) expectation->argStruct;
506 //newObj->data = oo->expectation->data;
508 oo->computeFun = ba81Compute;
509 oo->setVarGroup = ba81SetFreeVarGroup;
510 oo->destructFun = ba81Destroy;
511 oo->gradientAvailable = TRUE;
512 oo->hessianAvailable = TRUE;
515 PROTECT(tmp = GET_SLOT(rObj, install("rescale")));
516 state->rescale = asLogical(tmp);
519 omxNewMatrixFromSlot(rObj, globalState, "ItemParam");
521 if (estate->EitemParam->rows != state->itemParam->rows ||
522 estate->EitemParam->cols != state->itemParam->cols) {
523 error("ItemParam and EItemParam matrices must be the same dimension");
527 omxNewMatrixFromSlot(rObj, globalState, "CustomPrior");
529 int maxParam = state->itemParam->rows;
530 state->derivPadSize = maxParam + maxParam*(1+maxParam)/2;
532 state->tmpLatentMean = Realloc(NULL, estate->maxDims, double);
533 state->tmpLatentCov = Realloc(NULL, estate->maxDims * estate->maxDims, double);
535 int numItems = state->itemParam->cols;
536 for (int ix=0; ix < numItems; ix++) {
537 double *spec = omxMatrixColumn(estate->itemSpec, ix);
538 int id = spec[RPF_ISpecID];
539 if (id < 0 || id >= rpf_numModels) {
540 error("ItemSpec column %d has unknown item model %d", ix, id);