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 "omxFitFunction.h"
  19 #include "omxExpectationBA81.h"
  20 #include "omxOpenmpWrap.h"
  21 #include "libifa-rpf.h"
  22
  23 static const char *NAME = "FitFunctionBA81";
  24
  25 struct BA81FitState {
  26
  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;
  34         bool rescale;
  35         omxMatrix *customPrior;
  36         int choleskyError;
  37         double *tmpLatentMean;    // maxDims
  38         double *tmpLatentCov;     // maxDims * maxDims ; only lower triangle is used
  39         int fitCount;
  40         int gradientCount;
  41
  42         std::vector< FreeVarGroup* > varGroups;
  43         FreeVarGroup *latentFVG;
  44
  45         BA81FitState();
  46 };
  47
  48 BA81FitState::BA81FitState()
  49 {
  50         itemParam = NULL;
  51         thrDeriv = NULL;
  52         paramMap = NULL;
  53         latentFVG = NULL;
  54         customPrior = NULL;
  55         fitCount = 0;
  56         gradientCount = 0;
  57         tmpLatentMean = NULL;
  58         tmpLatentCov = NULL;
  59 }
  60
  61 static void buildParamMap(omxFitFunction* oo)
  62 {
  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;
  71
  72         state->latentMeanMap.assign(maxAbilities, -1);
  73         state->latentCovMap.assign(maxAbilities * maxAbilities, -1);
  74
  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;
  84                         }
  85                 }
  86         }
  87
  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;
  91         }
  92
  93         size_t numFreeParams = oo->freeVarGroup->vars.size();
  94         int *pRow = Realloc(NULL, numFreeParams, int);
  95         int *pCol = Realloc(NULL, numFreeParams, int);
  96
  97         for (size_t px=0; px < numFreeParams; px++) {
  98                 pRow[px] = -1;
  99                 pCol[px] = -1;
 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) {
 105                                 pRow[px] = loc->row;
 106                                 pCol[px] = loc->col;
 107                                 int at = pCol[px] * state->derivPadSize + pRow[px];
 108                                 state->paramMap[at] = px;
 109
 110                                 const double *spec = omxMatrixColumn(estate->itemSpec, loc->col);
 111                                 int id = spec[RPF_ISpecID];
 112                                 double upper, lower;
 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;
 116                         }
 117                 }
 118         }
 119
 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);
 126                         int r1 = pRow[p1];
 127                         int r2 = pRow[p2];
 128                         if (r1 > r2) { int tmp=r1; r1=r2; r2=tmp; }
 129                         int rowOffset = 0;
 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);
 134                 }
 135         }
 136
 137         Free(pRow);
 138         Free(pCol);
 139
 140         state->thrDeriv = Realloc(NULL, itemParam->cols * state->derivPadSize * Global->numThreads, double);
 141 }
 142
 143 OMXINLINE static double
 144 ba81Fit1Ordinate(omxFitFunction* oo, const int *quad, const double *weight, int want)
 145 {
 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);
 155
 156         double where[maxDims];
 157         pointToWhere(estate, quad, where, maxDims);
 158
 159         double *outcomeProb = computeRPF(estate, itemParam, quad); // avoid malloc/free? TODO
 160         if (!outcomeProb) return 0;
 161
 162         double thr_ll = 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];
 167
 168                 double area = exp(logAreaProduct(estate, quad, estate->Sgroup[ix]));   // avoid exp() here? TODO
 169                 for (int ox=0; ox < iOutcomes; ox++) {
 170 #if 0
 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");
 176                         }
 177 #endif
 178                         double got = weight[ox] * outcomeProb[ix * maxOutcomes + ox];
 179                         thr_ll += got * area;
 180                 }
 181
 182                 if (do_deriv) {
 183                         double *iparam = omxMatrixColumn(itemParam, ix);
 184                         double *pad = myDeriv + ix * state->derivPadSize;
 185                         (*rpf_model[id].dLL1)(spec, iparam, where, area, weight, pad);
 186                 }
 187                 weight += iOutcomes;
 188         }
 189
 190         Free(outcomeProb);
 191
 192         return thr_ll;
 193 }
 194
 195 static double
 196 ba81ComputeMFit1(omxFitFunction* oo, int want, double *gradient, double *hessian)
 197 {
 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;
 205
 206         double ll = 0;
 207         if (customPrior) {
 208                 omxRecompute(customPrior);
 209                 ll = customPrior->data[0];
 210                 // need deriv adjustment TODO
 211         }
 212
 213         if (!isfinite(ll)) {
 214                 omxPrint(itemParam, "item param");
 215                 error("Bayesian prior returned %g; do you need to add a lbound/ubound?", ll);
 216         }
 217
 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
 221                 int quad[maxDims];
 222                 decodeLocation(qx, maxDims, estate->quadGridSize, quad);
 223                 double *weight = estate->expected + qx * totalOutcomes;
 224                 double thr_ll = ba81Fit1Ordinate(oo, quad, weight, want);
 225
 226 #pragma omp atomic
 227                 ll += thr_ll;
 228         }
 229
 230         if (gradient) {
 231                 double *deriv0 = state->thrDeriv;
 232
 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];
 237                 }
 238
 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);
 246                 }
 247
 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;
 253
 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]);
 267                         }
 268
 269                         if (to < numFreeParams) {
 270                                 gradient[to] -= deriv0[ox];
 271                         } else {
 272                                 hessian[to - numFreeParams] -= deriv0[ox];
 273                         }
 274                 }
 275         }
 276
 277         return -ll;
 278 }
 279
 280 static void
 281 moveLatentDistribution(omxFitFunction *oo, FitContext *fc,
 282                        double *ElatentMean, double *ElatentCov)
 283 {
 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;
 293
 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];
 300                 int idx = 0;
 301                 for (int dx=0; dx < design->rows; dx++) {
 302                         if (isfinite(rawDesign[dx])) {
 303                                 idesign[idx++] = rawDesign[dx]-1;
 304                         } else {
 305                                 idesign[idx++] = -1;
 306                         }
 307                 }
 308                 for (int d1=0; d1 < idx; d1++) {
 309                         if (idesign[d1] == -1) {
 310                                 tmpLatentMean[d1] = 0;
 311                         } else {
 312                                 tmpLatentMean[d1] = ElatentMean[idesign[d1]];
 313                         }
 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;
 318                                 } else {
 319                                         tmpLatentCov[d2 * maxDims + d1] = ElatentCov[cell];
 320                                 }
 321                         }
 322                 }
 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("");
 329                                 }
 330                         }
 331                 }
 332                 double *iparam = omxMatrixColumn(itemParam, ix);
 333                 int *mask = state->paramMap + state->derivPadSize * ix;
 334                 rpf_model[id].rescale(spec, iparam, mask, tmpLatentMean, tmpLatentCov);
 335         }
 336
 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);
 343                         }
 344                 }
 345         }
 346 }
 347
 348 static void
 349 schilling_bock_2005_rescale(omxFitFunction *oo, FitContext *fc)
 350 {
 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;
 356
 357         //mxLog("schilling bock\n");
 358         //pda(ElatentMean, maxAbilities, 1);
 359         //pda(ElatentCov, maxAbilities, maxAbilities);
 360         //omxPrint(design, "design");
 361
 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?
 367                 return;
 368         }
 369
 370         moveLatentDistribution(oo, fc, ElatentMean, ElatentCov);
 371         fc->copyParamToModel(globalState);
 372 }
 373
 374 OMXINLINE static void
 375 updateLatentParam(omxFitFunction* oo, FitContext *fc)
 376 {
 377         BA81FitState *state = (BA81FitState*) oo->argStruct;
 378         BA81Expect *estate = (BA81Expect*) oo->expectation->argStruct;
 379         int maxAbilities = estate->maxAbilities;
 380
 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;
 389                         fc->est[vx] = val;
 390                 }
 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;
 399                         fc->est[vx] = val;
 400                 }
 401         }
 402
 403         fc->copyParamToModel(globalState);
 404 }
 405
 406 void ba81SetFreeVarGroup(omxFitFunction *oo, FreeVarGroup *fvg) // too ad hoc? TODO
 407 {
 408         if (!oo->argStruct) { // ugh!
 409                 BA81FitState *state = new BA81FitState;
 410                 oo->argStruct = state;
 411         }
 412
 413         BA81FitState *state = (BA81FitState*) oo->argStruct;
 414
 415         state->varGroups.push_back(fvg);
 416         if (state->varGroups.size() == 2) {
 417                 int small = 0;
 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())
 421                         small = 1;
 422                 oo->freeVarGroup = state->varGroups[small];
 423                 state->latentFVG = state->varGroups[!small];
 424         } else if (state->varGroups.size() > 2) {
 425                 // ignore
 426         }
 427 }
 428
 429 static double
 430 ba81ComputeFit(omxFitFunction* oo, int want, FitContext *fc)
 431 {
 432         BA81FitState *state = (BA81FitState*) oo->argStruct;
 433
 434         if (!state->paramMap) buildParamMap(oo);
 435
 436         if (want & FF_COMPUTE_PREOPTIMIZE) {
 437                 if (state->rescale) schilling_bock_2005_rescale(oo, fc); // how does this work in multigroup? TODO
 438                 return 0;
 439         }
 440
 441         BA81Expect *estate = (BA81Expect*) oo->expectation->argStruct;
 442
 443         ++state->fitCount;
 444
 445         if (estate->type == EXPECTATION_AUGMENTED) {
 446                 if (fc->varGroup != oo->freeVarGroup) error("FreeVarGroup mismatch");
 447
 448                 if (want & FF_COMPUTE_GRADIENT) ++state->gradientCount;
 449
 450                 omxMatrix *itemParam = state->itemParam;
 451                 OMXZERO(state->thrDeriv, state->derivPadSize * itemParam->cols * Global->numThreads);
 452
 453                 for (size_t nx=0; nx < state->NAtriangle.size(); ++nx) {
 454                         fc->hess[ state->NAtriangle[nx] ] = nan("symmetric");
 455                 }
 456
 457                 double got = ba81ComputeMFit1(oo, want, fc->grad, fc->hess);
 458                 return got;
 459         } else if (estate->type == EXPECTATION_OBSERVED) {
 460                 if (fc->varGroup != state->latentFVG) error("FreeVarGroup mismatch");
 461
 462                 updateLatentParam(oo, fc);
 463
 464                 double *patternLik = estate->patternLik;
 465                 int *numIdentical = estate->numIdentical;
 466                 int numUnique = estate->numUnique;
 467                 double got = 0;
 468                 for (int ux=0; ux < numUnique; ux++) {
 469                         got += numIdentical[ux] * patternLik[ux];
 470                 }
 471                 return -2 * got;
 472         }
 473 }
 474
 475 static void ba81Compute(omxFitFunction *oo, int want, FitContext *fc)
 476 {
 477         if (!want) return;
 478         oo->matrix->data[0] = ba81ComputeFit(oo, want, fc);
 479 }
 480
 481 static void ba81Destroy(omxFitFunction *oo) {
 482         BA81FitState *state = (BA81FitState *) oo->argStruct;
 483
 484         omxFreeAllMatrixData(state->customPrior);
 485         Free(state->paramMap);
 486         Free(state->thrDeriv);
 487         Free(state->tmpLatentMean);
 488         Free(state->tmpLatentCov);
 489         omxFreeAllMatrixData(state->itemParam);
 490         delete state;
 491 }
 492
 493 void omxInitFitFunctionBA81(omxFitFunction* oo)
 494 {
 495         if (!oo->argStruct) { // ugh!
 496                 BA81FitState *state = new BA81FitState;
 497                 oo->argStruct = state;
 498         }
 499
 500         BA81FitState *state = (BA81FitState*) oo->argStruct;
 501         SEXP rObj = oo->rObj;
 502
 503         omxExpectation *expectation = oo->expectation;
 504         BA81Expect *estate = (BA81Expect*) expectation->argStruct;
 505
 506         //newObj->data = oo->expectation->data;
 507
 508         oo->computeFun = ba81Compute;
 509         oo->setVarGroup = ba81SetFreeVarGroup;
 510         oo->destructFun = ba81Destroy;
 511         oo->gradientAvailable = TRUE;
 512         oo->hessianAvailable = TRUE;
 513
 514         SEXP tmp;
 515         PROTECT(tmp = GET_SLOT(rObj, install("rescale")));
 516         state->rescale = asLogical(tmp);
 517
 518         state->itemParam =
 519                 omxNewMatrixFromSlot(rObj, globalState, "ItemParam");
 520
 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");
 524         }
 525
 526         state->customPrior =
 527                 omxNewMatrixFromSlot(rObj, globalState, "CustomPrior");
 528
 529         int maxParam = state->itemParam->rows;
 530         state->derivPadSize = maxParam + maxParam*(1+maxParam)/2;
 531
 532         state->tmpLatentMean = Realloc(NULL, estate->maxDims, double);
 533         state->tmpLatentCov = Realloc(NULL, estate->maxDims * estate->maxDims, double);
 534
 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);
 541                 }
 542         }
 543 }