Fixed a bug where could not be accessed from a model object, and adjusted one test...
[openmx:openmx.git] / src / omxExpectationBA81.h
1 /*
2   Copyright 2012-2014 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 #include "ba81quad.h"
24
25 enum expectation_type {
26         EXPECTATION_AUGMENTED, // E-M
27         EXPECTATION_OBSERVED,  // regular
28 };
29
30 template <typename CovType>
31 struct BA81EstepBase {
32         void addRow1(class ifaGroup *state, int px, double *Qweight, double *out);
33 };
34
35 template <typename T, typename CovType>
36 struct BA81Estep : BA81EstepBase<CovType> {
37         std::vector<double> thrExpected;
38
39         void begin(class ifaGroup *state, T extraData);
40         void addRow(class ifaGroup *state, T extraData, int px, double *Qweight, int thrId);
41         void recordTable(class ifaGroup *state, T extraData);
42         bool hasEnd() { return true; }
43 };
44
45 template <typename T>
46 struct BA81LatentFixed {
47         void begin(class ifaGroup *state, T extraData) {}
48         void normalizeWeights(class ifaGroup *state, T extraData, int px, double *Qweight, double weight, int thrid);
49         void end(class ifaGroup *state, T extraData) {};
50         bool hasEnd() { return false; }
51 };
52
53 template <typename T>
54 struct BA81LatentSummary {
55         void begin(class ifaGroup *state, T extraData);
56         void normalizeWeights(class ifaGroup *state, T extraData, int px, double *Qweight, double weight, int thrId);
57         void end(class ifaGroup *state, T extraData);
58         bool hasEnd() { return true; }
59 };
60
61 class BA81Expect {
62  public:
63         class ifaGroup grp;
64         int totalOutcomes() { return grp.totalOutcomes; }
65         const double *itemSpec(int ix) { return grp.spec[ix]; }
66         int numItems() { return grp.numItems(); }
67         int getNumUnique() { return (int) grp.rowMap.size(); }
68         int itemOutcomes(int ix) { return grp.itemOutcomes[ix]; }
69
70         const char *name;              // from omxExpectation
71         double LogLargestDouble;       // should be const but need constexpr
72         double LargestDouble;          // should be const but need constexpr
73
74         // data characteristics
75         omxData *data;
76         double weightSum;                // sum of rowWeight
77         // aggregate distribution of data in quadrature
78         std::vector<double> thrDweight;  // quad.weightTableSize * numThreads
79
80         // quadrature related
81         struct ba81NormalQuad &getQuad() { return grp.quad; }
82
83         // estimation related
84         omxMatrix *itemParam;
85         double *EitemParam;
86         double SmallestPatternLik;
87         double *expected;                     // totalOutcomes * totalQuadPoints (E-step table)
88         bool expectedUsed;
89         int ElatentVersion;
90
91         omxMatrix *_latentMeanOut;
92         omxMatrix *_latentCovOut;
93         template <typename Tmean, typename Tcov>
94         void getLatentDistribution(FitContext *fc, Eigen::MatrixBase<Tmean> &mean, Eigen::MatrixBase<Tcov> &cov);
95
96         omxMatrix *estLatentMean;
97         omxMatrix *estLatentCov;
98
99         int itemParamVersion;
100         int latentParamVersion;
101         enum expectation_type type;
102         int verbose;
103         bool debugInternal;
104         struct omxFitFunction *fit;  // weak pointer
105
106         BA81Expect() : grp(Global->numThreads, true) {};
107         const char *getLatentIncompatible(BA81Expect *other);
108 };
109
110 template <typename Tmean, typename Tcov>
111 void BA81Expect::getLatentDistribution(FitContext *fc, Eigen::MatrixBase<Tmean> &mean, Eigen::MatrixBase<Tcov> &cov)
112 {
113         mean.derived().resize(grp.maxAbilities);
114         if (!_latentMeanOut) {
115                 mean.setZero();
116         } else {
117                 omxRecompute(_latentMeanOut, fc);
118                 memcpy(mean.derived().data(), _latentMeanOut->data, sizeof(double) * grp.maxAbilities);
119         }
120         
121         cov.derived().resize(grp.maxAbilities, grp.maxAbilities);
122         if (!_latentCovOut) {
123                 cov.setIdentity();
124         } else {
125                 omxRecompute(_latentCovOut, fc);
126                 memcpy(cov.derived().data(), _latentCovOut->data, sizeof(double) * grp.maxAbilities * grp.maxAbilities);
127         }
128 }
129
130 extern const struct rpf *rpf_model;
131 extern int rpf_numModels;
132
133 OMXINLINE static void
134 gramProduct(double *vec, size_t len, double *out)
135 {
136         int cell = 0;
137         for (size_t v1=0; v1 < len; ++v1) {
138                 for (size_t v2=0; v2 <= v1; ++v2) {
139                         out[cell] = vec[v1] * vec[v2];
140                         ++cell;
141                 }
142         }
143 }
144
145 void ba81SetupQuadrature(omxExpectation* oo);
146
147 void ba81AggregateDistributions(std::vector<struct omxExpectation *> &expectation,
148                                 int *version, omxMatrix *meanMat, omxMatrix *covMat);
149
150 static const double BA81_MIN_VARIANCE = .01;
151
152 #endif