Reduce rpf evaluations in two-tier models
[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 enum expectation_type {
30         EXPECTATION_UNINITIALIZED,
31         EXPECTATION_AUGMENTED, // E-M
32         EXPECTATION_OBSERVED,  // regular
33 };
34
35 typedef struct {
36
37         // data characteristics
38         omxData *data;
39         int numUnique;
40         int *numIdentical;        // length numUnique
41         double *logNumIdentical;  // length numUnique
42         int *rowMap;              // length numUnique, index of first instance of pattern
43
44         // item description related
45         omxMatrix *itemSpec;
46         int maxOutcomes;
47         int maxDims;
48         int maxAbilities;
49         int numSpecific;
50         int *Sgroup;              // item's specific group 0..numSpecific-1
51         omxMatrix *design;        // items * maxDims
52
53         // quadrature related
54         double Qwidth;
55         double targetQpoints;
56         long quadGridSize;
57         long totalQuadPoints;     // quadGridSize ^ maxDims
58         long totalPrimaryPoints;  // totalQuadPoints except for specific dim TODO
59         std::vector<double> Qpoint;           // quadGridSize
60         std::vector<double> priLogQarea;      // totalPrimaryPoints
61         std::vector<double> speLogQarea;      // quadGridSize * numSpecific
62
63         // estimation related
64         omxMatrix *EitemParam;    // E step version
65         int cacheLXK;
66         double *lxk;              // wo/cache, numUnique * thread
67         std::vector<double> allElxk;          // numUnique * thread
68         std::vector<double> Eslxk;            // numUnique * #specific dimensions * thread
69         double *patternLik;       // numUnique
70         int totalOutcomes;
71         double *expected;         // totalOutcomes * totalQuadPoints
72         double *ElatentMean;      // maxAbilities * numUnique
73         double *ElatentCov;       // maxAbilities * maxAbilities * numUnique ; only lower triangle is used
74         omxMatrix *latentMeanOut;
75         omxMatrix *latentCovOut;
76
77         enum expectation_type type;
78         enum score_option scores;
79 } BA81Expect;
80
81 extern const struct rpf *rpf_model;
82 extern int rpf_numModels;
83
84 double *computeRPF(BA81Expect *state, omxMatrix *itemParam, const int *quad);
85
86 OMXINLINE static void
87 pointToWhere(BA81Expect *state, const int *quad, double *where, int upto)
88 {
89         for (int dx=0; dx < upto; dx++) {
90                 where[dx] = state->Qpoint[quad[dx]];
91         }
92 }
93
94 OMXINLINE static long
95 encodeLocation(const int dims, const long grid, const int *quad)
96 {
97         long qx = 0;
98         for (int dx=dims-1; dx >= 0; dx--) {
99                 qx = qx * grid;
100                 qx += quad[dx];
101         }
102         return qx;
103 }
104
105 OMXINLINE static void
106 decodeLocation(long qx, const int dims, const long grid, int *quad)
107 {
108         for (int dx=0; dx < dims; dx++) {
109                 quad[dx] = qx % grid;
110                 qx = qx / grid;
111         }
112 }
113
114 OMXINLINE static double
115 logAreaProduct(BA81Expect *state, const int *quad, const int sg)
116 {
117         int maxDims = state->maxDims;
118         if (state->numSpecific == 0) {
119                 long qloc = encodeLocation(maxDims, state->quadGridSize, quad);
120                 return state->priLogQarea[qloc];
121         } else {
122                 long priloc = encodeLocation(maxDims-1, state->quadGridSize, quad);
123                 return (state->priLogQarea[priloc] +
124                         state->speLogQarea[sg * state->quadGridSize + quad[maxDims - 1]]);
125         }
126 }
127
128 // debug tools
129 void pda(const double *ar, int rows, int cols);
130 void pia(const int *ar, int rows, int cols);
131
132 #endif