Move EAP scores to expectation slot scores.out
[openmx:openmx.git] / src / omxExpectationBA81.h
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 #ifndef _OMX_EXPECTATIONBA81_H_
19 #define _OMX_EXPECTATIONBA81_H_
20
21 #include "omxExpectation.h"
22
23 enum score_option {
24         SCORES_OMIT,
25         SCORES_UNIQUE,
26         SCORES_FULL
27 };
28
29 typedef struct {
30
31         // data characteristics
32         omxData *data;
33         int numUnique;
34         int *numIdentical;        // length numUnique
35         double *logNumIdentical;  // length numUnique
36         int *rowMap;              // length numUnique
37
38         // item description related
39         omxMatrix *itemSpec;
40         int maxOutcomes;
41         int maxDims;
42         int maxAbilities;
43         int numSpecific;
44         int *Sgroup;              // item's specific group 0..numSpecific-1
45         omxMatrix *design;        // items * maxDims
46
47         // quadrature related
48         double Qwidth;
49         double targetQpoints;
50         long quadGridSize;
51         long totalQuadPoints;     // quadGridSize ^ maxDims
52         long totalPrimaryPoints;  // totalQuadPoints except for specific dim TODO
53         double *Qpoint;           // quadGridSize
54         double *priLogQarea;      // totalPrimaryPoints
55         double *speLogQarea;      // quadGridSize * numSpecific
56
57         // estimation related
58         omxMatrix *EitemParam;    // E step version
59         int cacheLXK;
60         double *lxk;              // wo/cache, numUnique * thread
61         double *allSlxk;          // numUnique * thread
62         double *Slxk;             // numUnique * #specific dimensions * thread
63         double *patternLik;       // numUnique
64         int totalOutcomes;
65         double *expected;         // totalOutcomes * totalQuadPoints
66         double *ElatentMean;      // maxAbilities * numUnique
67         double *ElatentCov;       // maxAbilities * maxAbilities * numUnique ; only lower triangle is used
68         omxMatrix *latentMeanOut;
69         omxMatrix *latentCovOut;
70
71         enum score_option scores;
72 } BA81Expect;
73
74 extern const struct rpf *rpf_model;
75 extern int rpf_numModels;
76
77 double *
78 computeRPF(omxMatrix *itemSpec, omxMatrix *design, omxMatrix *itemParam,
79            int maxDims, int maxOutcomes, const int *quad, const double *Qpoint);
80
81 OMXINLINE static void
82 pointToWhere(const double *Qpoint, const int *quad, double *where, int upto)
83 {
84         for (int dx=0; dx < upto; dx++) {
85                 where[dx] = Qpoint[quad[dx]];
86         }
87 }
88
89 OMXINLINE static long
90 encodeLocation(const int dims, const long grid, const int *quad)
91 {
92         long qx = 0;
93         for (int dx=dims-1; dx >= 0; dx--) {
94                 qx = qx * grid;
95                 qx += quad[dx];
96         }
97         return qx;
98 }
99
100 OMXINLINE static void
101 decodeLocation(long qx, const int dims, const long grid, int *quad)
102 {
103         for (int dx=0; dx < dims; dx++) {
104                 quad[dx] = qx % grid;
105                 qx = qx / grid;
106         }
107 }
108
109 OMXINLINE static double
110 logAreaProduct(BA81Expect *state, const int *quad, const int sg)
111 {
112         int maxDims = state->maxDims;
113         if (state->numSpecific == 0) {
114                 long qloc = encodeLocation(maxDims, state->quadGridSize, quad);
115                 return state->priLogQarea[qloc];
116         } else {
117                 long priloc = encodeLocation(maxDims-1, state->quadGridSize, quad);
118                 return (state->priLogQarea[priloc] +
119                         state->speLogQarea[sg * state->quadGridSize + quad[maxDims - 1]]);
120         }
121 }
122
123 // debug tools
124 void pda(const double *ar, int rows, int cols);
125 void pia(const int *ar, int rows, int cols);
126
127 #endif