Transpose outcomeProb for better cache behavior
[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 #include "omxOpenmpWrap.h"
23
24 enum score_option {
25         SCORES_OMIT,
26         SCORES_UNIQUE,
27         SCORES_FULL
28 };
29
30 enum expectation_type {
31         EXPECTATION_UNINITIALIZED,
32         EXPECTATION_AUGMENTED, // E-M
33         EXPECTATION_OBSERVED,  // regular
34 };
35
36 struct BA81Expect {
37         double LogLargestDouble;       // should be const but need constexpr
38         double LargestDouble;          // should be const but need constexpr
39         double OneOverLargestDouble;   // should be const but need constexpr
40
41         // data characteristics
42         omxData *data;
43         int numUnique;
44         int *numIdentical;        // length numUnique
45         int *rowMap;              // length numUnique, index of first instance of pattern
46
47         // item description related
48         std::vector<const double*> itemSpec;
49         std::vector<int> itemOutcomes;
50         int maxDims;
51         int maxAbilities;
52         int numSpecific;
53         int *Sgroup;              // item's specific group 0..numSpecific-1
54         omxMatrix *design;        // items * maxDims
55
56         // quadrature related
57         double Qwidth;
58         double targetQpoints;
59         long quadGridSize;
60         long totalQuadPoints;     // quadGridSize ^ maxDims
61         long totalPrimaryPoints;  // totalQuadPoints except for specific dim TODO
62         std::vector<double> Qpoint;           // quadGridSize
63         std::vector<double> priQarea;         // totalPrimaryPoints
64         std::vector<double> speQarea;         // quadGridSize * numSpecific
65
66         // estimation related
67         omxMatrix *customPrior;
68         omxMatrix *itemParam;
69         int cacheLXK;
70         bool LXKcached;
71         double *lxk;              // wo/cache, numUnique * thread
72         double *allElxk;          // numUnique * thread
73         double *Eslxk;            // numUnique * #specific dimensions * thread
74         double *patternLik;       // numUnique
75         double SmallestPatternLik;
76         int excludedPatterns;
77         int totalOutcomes;
78         double *outcomeProb;      // totalOutcomes * totalQuadPoints
79         double *expected;         // totalOutcomes * totalQuadPoints
80         std::vector<double> ElatentMean;      // maxAbilities
81         std::vector<double> ElatentCov;       // maxAbilities * maxAbilities
82         omxMatrix *latentMeanOut;
83         omxMatrix *latentCovOut;
84         double *EiCache;          // totalPrimaryPoints * numUnique (two-tier only)
85
86         int itemParamVersion;
87         int latentParamVersion;
88         enum expectation_type type;
89         enum score_option scores;
90         bool verbose;
91 };
92
93 extern const struct rpf *rpf_model;
94 extern int rpf_numModels;
95
96 template<typename _Tp> class omxBuffer {
97         typedef _Tp value_type;
98         typedef _Tp *pointer;
99         typedef _Tp &reference;
100         typedef const _Tp &const_reference;
101         pointer _M_start;
102  public:
103         omxBuffer(size_t size) {
104                 _M_start = new _Tp[size];
105         }
106         omxBuffer(size_t size, const value_type &value) {
107                 _M_start = new _Tp[size];
108                 std::fill_n(_M_start, size, value);
109         }
110         ~omxBuffer() {
111                 delete [] _M_start;
112         }
113         reference operator[](size_t __n)
114         { return *(_M_start + __n); }
115         const_reference operator[](size_t __n) const
116         { return *(_M_start + __n); }
117         pointer data() const
118         { return _M_start; }
119 };
120
121 OMXINLINE static int
122 triangleLoc1(int diag)
123 {
124         //if (diag < 1) error("Out of domain");
125         return (diag) * (diag+1) / 2;   // 0 1 3 6 10 15 ..
126 }
127
128 OMXINLINE static int
129 triangleLoc0(int diag)
130 {
131         //if (diag < 0) error("Out of domain");
132         return triangleLoc1(diag+1) - 1;  // 0 2 5 9 14 ..
133 }
134
135 OMXINLINE static double *
136 eBase(BA81Expect *state, int thr)
137 {
138         return state->allElxk + thr * state->numUnique;
139 }
140
141 OMXINLINE static double *
142 esBase(BA81Expect *state, int thr)
143 {
144         return state->Eslxk + thr * state->numUnique * state->numSpecific;
145 }
146
147 OMXINLINE static void
148 pointToWhere(BA81Expect *state, const int *quad, double *where, int upto)
149 {
150         for (int dx=0; dx < upto; dx++) {
151                 where[dx] = state->Qpoint[quad[dx]];
152         }
153 }
154
155 OMXINLINE static void
156 decodeLocation(long qx, const int dims, const long grid, int *quad)
157 {
158         for (int dx=dims-1; dx >= 0; --dx) {
159                 quad[dx] = qx % grid;
160                 qx = qx / grid;
161         }
162 }
163
164 OMXINLINE static double
165 areaProduct(BA81Expect *state, long qx, int sx, const int sg)
166 {
167         if (state->numSpecific == 0) {
168                 return state->priQarea[qx];
169         } else {
170                 if (sx == -1) {
171                         sx = qx % state->quadGridSize;
172                         qx /= state->quadGridSize;
173                 }
174                 return (state->priQarea[qx] * state->speQarea[sg * state->quadGridSize + sx]);
175         }
176 }
177
178 OMXINLINE static void
179 gramProduct(double *vec, size_t len, double *out)
180 {
181         int cell = 0;
182         for (size_t v1=0; v1 < len; ++v1) {
183                 for (size_t v2=0; v2 <= v1; ++v2) {
184                         out[cell] = vec[v1] * vec[v2];
185                         ++cell;
186                 }
187         }
188 }
189
190 OMXINLINE static bool
191 validPatternLik(BA81Expect *state, double pl)
192 {
193         return isfinite(pl) && pl > state->SmallestPatternLik;
194 }
195
196 // debug tools
197 void pda(const double *ar, int rows, int cols);
198 void pia(const int *ar, int rows, int cols);
199
200 #endif