84 | | class RV : public root { |
85 | | protected: |
86 | | //! size of the data vector |
87 | | int dsize; |
88 | | //! number of individual rvs |
89 | | int len; |
90 | | //! Vector of unique IDs |
91 | | ivec ids; |
92 | | //! Vector of shifts from current time |
93 | | ivec times; |
94 | | |
95 | | private: |
96 | | //! auxiliary function used in constructor |
97 | | void init ( const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times ); |
98 | | int init ( const string &name, int size ); |
99 | | public: |
100 | | //! \name Constructors |
101 | | //!@{ |
102 | | |
103 | | //! Full constructor |
104 | | RV ( const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times ) { |
105 | | init ( in_names, in_sizes, in_times ); |
106 | | } |
107 | | |
108 | | //! Constructor with times=0 |
109 | | RV ( const Array<std::string> &in_names, const ivec &in_sizes ) { |
110 | | init ( in_names, in_sizes, zeros_i ( in_names.length() ) ); |
111 | | } |
112 | | |
113 | | //! Constructor with sizes=1, times=0 |
114 | | RV ( const Array<std::string> &in_names ) { |
115 | | init ( in_names, ones_i ( in_names.length() ), zeros_i ( in_names.length() ) ); |
116 | | } |
117 | | |
118 | | //! Constructor of empty RV |
119 | | RV() : dsize ( 0 ), len ( 0 ), ids ( 0 ), times ( 0 ) {} |
120 | | //! Constructor of a single RV with given id |
121 | | RV ( string name, int sz, int tm = 0 ); |
122 | | //!@} |
123 | | |
124 | | //! \name Access functions |
125 | | //!@{ |
126 | | |
127 | | //! State output, e.g. for debugging. |
128 | | friend std::ostream &operator<< ( std::ostream &os, const RV &rv ); |
129 | | |
130 | | int _dsize() const { |
131 | | return dsize; |
132 | | } |
133 | | |
134 | | //! Recount size of the corresponding data vector |
135 | | int countsize() const; |
136 | | ivec cumsizes() const; |
137 | | int length() const { |
138 | | return len; |
139 | | } |
140 | | int id ( int at ) const { |
141 | | return ids ( at ); |
142 | | } |
143 | | int size ( int at ) const { |
144 | | return RV_SIZES ( ids ( at ) ); |
145 | | } |
146 | | int time ( int at ) const { |
147 | | return times ( at ); |
148 | | } |
149 | | std::string name ( int at ) const { |
150 | | return RV_NAMES ( ids ( at ) ); |
151 | | } |
152 | | void set_time ( int at, int time0 ) { |
153 | | times ( at ) = time0; |
154 | | } |
155 | | //!@} |
156 | | |
157 | | //TODO why not inline and later?? |
158 | | |
159 | | //! \name Algebra on Random Variables |
160 | | //!@{ |
161 | | |
162 | | //! Find indices of self in another rv, \return ivec of the same size as self. |
163 | | ivec findself ( const RV &rv2 ) const; |
164 | | //! Compare if \c rv2 is identical to this \c RV |
165 | | bool equal ( const RV &rv2 ) const; |
166 | | //! Add (concat) another variable to the current one, \return true if all rv2 were added, false if rv2 is in conflict |
167 | | bool add ( const RV &rv2 ); |
168 | | //! Subtract another variable from the current one |
169 | | RV subt ( const RV &rv2 ) const; |
170 | | //! Select only variables at indices ind |
171 | | RV subselect ( const ivec &ind ) const; |
172 | | |
173 | | //! Select only variables at indices ind |
174 | | RV operator() ( const ivec &ind ) const { |
175 | | return subselect ( ind ); |
176 | | } |
177 | | |
178 | | //! Select from data vector starting at di1 to di2 |
179 | | RV operator() ( int di1, int di2 ) const; |
180 | | |
181 | | //! Shift \c time by delta. |
182 | | void t ( int delta ); |
183 | | //!@} |
184 | | |
185 | | //!\name Relation to vectors |
186 | | //!@{ |
187 | | |
188 | | //! generate \c str from rv, by expanding sizes TODO to_string.. |
189 | | str tostr() const; |
190 | | //! when this rv is a part of bigger rv, this function returns indices of self in the data vector of the bigger crv. |
191 | | //! Then, data can be copied via: data_of_this = cdata(ind); |
192 | | ivec dataind ( const RV &crv ) const; |
193 | | //! generate mutual indices when copying data between self and crv. |
194 | | //! Data are copied via: data_of_this(selfi) = data_of_rv2(rv2i) |
195 | | void dataind ( const RV &rv2, ivec &selfi, ivec &rv2i ) const; |
196 | | //! Minimum time-offset |
197 | | int mint() const { |
198 | | return min ( times ); |
199 | | } |
200 | | //!@} |
201 | | |
202 | | // TODO aktualizovat dle soucasneho UI |
203 | | /*! \brief UI for class RV (description of data vectors) |
204 | | |
205 | | \code |
206 | | rv = { |
207 | | type = "rv"; //identifier of the description |
208 | | // UNIQUE IDENTIFIER same names = same variable |
209 | | names = ["a", "b", "c", ...]; // which will be used e.g. in loggers |
210 | | |
211 | | //optional arguments |
212 | | sizes = [1, 2, 3, ...]; // (optional) default = ones() |
213 | | times = [-1, -2, 0, ...]; // time shifts with respect to current time (optional) default = zeros() |
214 | | } |
215 | | \endcode |
216 | | */ |
217 | | void from_setting ( const Setting &set ); |
218 | | |
219 | | // TODO dodelat void to_setting( Setting &set ) const; |
220 | | |
221 | | //! Invalidate all named RVs. Use before initializing any RV instances, with care... |
222 | | static void clear_all(); |
223 | | }; |
224 | | UIREGISTER ( RV ); |
| 86 | class RV : public root |
| 87 | { |
| 88 | protected: |
| 89 | //! size of the data vector |
| 90 | int dsize; |
| 91 | //! number of individual rvs |
| 92 | int len; |
| 93 | //! Vector of unique IDs |
| 94 | ivec ids; |
| 95 | //! Vector of shifts from current time |
| 96 | ivec times; |
| 97 | |
| 98 | private: |
| 99 | //! auxiliary function used in constructor |
| 100 | void init (const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times); |
| 101 | int init (const string &name, int size); |
| 102 | public: |
| 103 | //! \name Constructors |
| 104 | //!@{ |
| 105 | |
| 106 | //! Full constructor |
| 107 | RV (const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times) { |
| 108 | init (in_names, in_sizes, in_times); |
| 109 | } |
| 110 | |
| 111 | //! Constructor with times=0 |
| 112 | RV (const Array<std::string> &in_names, const ivec &in_sizes) { |
| 113 | init (in_names, in_sizes, zeros_i (in_names.length())); |
| 114 | } |
| 115 | |
| 116 | //! Constructor with sizes=1, times=0 |
| 117 | RV (const Array<std::string> &in_names) { |
| 118 | init (in_names, ones_i (in_names.length()), zeros_i (in_names.length())); |
| 119 | } |
| 120 | |
| 121 | //! Constructor of empty RV |
| 122 | RV() : dsize (0), len (0), ids (0), times (0) {} |
| 123 | //! Constructor of a single RV with given id |
| 124 | RV (string name, int sz, int tm = 0); |
| 125 | //!@} |
| 126 | |
| 127 | //! \name Access functions |
| 128 | //!@{ |
| 129 | |
| 130 | //! State output, e.g. for debugging. |
| 131 | friend std::ostream &operator<< (std::ostream &os, const RV &rv); |
| 132 | |
| 133 | int _dsize() const { |
| 134 | return dsize; |
| 135 | } |
| 136 | |
| 137 | //! Recount size of the corresponding data vector |
| 138 | int countsize() const; |
| 139 | ivec cumsizes() const; |
| 140 | int length() const { |
| 141 | return len; |
| 142 | } |
| 143 | int id (int at) const { |
| 144 | return ids (at); |
| 145 | } |
| 146 | int size (int at) const { |
| 147 | return RV_SIZES (ids (at)); |
| 148 | } |
| 149 | int time (int at) const { |
| 150 | return times (at); |
| 151 | } |
| 152 | std::string name (int at) const { |
| 153 | return RV_NAMES (ids (at)); |
| 154 | } |
| 155 | void set_time (int at, int time0) { |
| 156 | times (at) = time0; |
| 157 | } |
| 158 | //!@} |
| 159 | |
| 160 | //TODO why not inline and later?? |
| 161 | |
| 162 | //! \name Algebra on Random Variables |
| 163 | //!@{ |
| 164 | |
| 165 | //! Find indices of self in another rv, \return ivec of the same size as self. |
| 166 | ivec findself (const RV &rv2) const; |
| 167 | //! Compare if \c rv2 is identical to this \c RV |
| 168 | bool equal (const RV &rv2) const; |
| 169 | //! Add (concat) another variable to the current one, \return true if all rv2 were added, false if rv2 is in conflict |
| 170 | bool add (const RV &rv2); |
| 171 | //! Subtract another variable from the current one |
| 172 | RV subt (const RV &rv2) const; |
| 173 | //! Select only variables at indices ind |
| 174 | RV subselect (const ivec &ind) const; |
| 175 | |
| 176 | //! Select only variables at indices ind |
| 177 | RV operator() (const ivec &ind) const { |
| 178 | return subselect (ind); |
| 179 | } |
| 180 | |
| 181 | //! Select from data vector starting at di1 to di2 |
| 182 | RV operator() (int di1, int di2) const; |
| 183 | |
| 184 | //! Shift \c time by delta. |
| 185 | void t (int delta); |
| 186 | //!@} |
| 187 | |
| 188 | //!\name Relation to vectors |
| 189 | //!@{ |
| 190 | |
| 191 | //! generate \c str from rv, by expanding sizes TODO to_string.. |
| 192 | str tostr() const; |
| 193 | //! when this rv is a part of bigger rv, this function returns indices of self in the data vector of the bigger crv. |
| 194 | //! Then, data can be copied via: data_of_this = cdata(ind); |
| 195 | ivec dataind (const RV &crv) const; |
| 196 | //! generate mutual indices when copying data between self and crv. |
| 197 | //! Data are copied via: data_of_this(selfi) = data_of_rv2(rv2i) |
| 198 | void dataind (const RV &rv2, ivec &selfi, ivec &rv2i) const; |
| 199 | //! Minimum time-offset |
| 200 | int mint() const { |
| 201 | return min (times); |
| 202 | } |
| 203 | //!@} |
| 204 | |
| 205 | // TODO aktualizovat dle soucasneho UI |
| 206 | /*! \brief UI for class RV (description of data vectors) |
| 207 | |
| 208 | \code |
| 209 | rv = { |
| 210 | type = "rv"; //identifier of the description |
| 211 | // UNIQUE IDENTIFIER same names = same variable |
| 212 | names = ["a", "b", "c", ...]; // which will be used e.g. in loggers |
| 213 | |
| 214 | //optional arguments |
| 215 | sizes = [1, 2, 3, ...]; // (optional) default = ones() |
| 216 | times = [-1, -2, 0, ...]; // time shifts with respect to current time (optional) default = zeros() |
| 217 | } |
| 218 | \endcode |
| 219 | */ |
| 220 | void from_setting (const Setting &set); |
| 221 | |
| 222 | // TODO dodelat void to_setting( Setting &set ) const; |
| 223 | |
| 224 | //! Invalidate all named RVs. Use before initializing any RV instances, with care... |
| 225 | static void clear_all(); |
| 226 | }; |
| 227 | UIREGISTER (RV); |
259 | | class epdf : public root { |
260 | | protected: |
261 | | //! dimension of the random variable |
262 | | int dim; |
263 | | //! Description of the random variable |
264 | | RV rv; |
265 | | |
266 | | public: |
267 | | /*! \name Constructors |
268 | | Construction of each epdf should support two types of constructors: |
269 | | \li empty constructor, |
270 | | \li copy constructor, |
271 | | |
272 | | The following constructors should be supported for convenience: |
273 | | \li constructor followed by calling \c set_parameters() |
274 | | \li constructor accepting random variables calling \c set_rv() |
275 | | |
276 | | All internal data structures are constructed as empty. Their values (including sizes) will be set by method \c set_parameters(). This way references can be initialized in constructors. |
277 | | @{*/ |
278 | | epdf() : dim ( 0 ), rv() {}; |
279 | | epdf ( const epdf &e ) : dim ( e.dim ), rv ( e.rv ) {}; |
280 | | epdf ( const RV &rv0 ) : dim ( rv0._dsize() ) { |
281 | | set_rv ( rv0 ); |
282 | | }; |
283 | | void set_parameters ( int dim0 ) { |
284 | | dim = dim0; |
285 | | } |
286 | | //!@} |
287 | | |
288 | | //! \name Matematical Operations |
289 | | //!@{ |
290 | | |
291 | | //! Returns a sample, \f$ x \f$ from density \f$ f_x()\f$ |
292 | | virtual vec sample() const { |
293 | | it_error ( "not implemneted" ); |
294 | | return vec ( 0 ); |
295 | | }; |
296 | | //! Returns N samples, \f$ [x_1 , x_2 , \ldots \ \f$ from density \f$ f_x(rv)\f$ |
297 | | virtual mat sample_m ( int N ) const; |
298 | | //! Compute log-probability of argument \c val |
299 | | //! In case the argument is out of suport return -Infinity |
300 | | virtual double evallog ( const vec &val ) const { |
301 | | it_error ( "not implemneted" ); |
302 | | return 0.0; |
303 | | }; |
304 | | //! Compute log-probability of multiple values argument \c val |
305 | | virtual vec evallog_m ( const mat &Val ) const { |
306 | | vec x ( Val.cols() ); |
307 | | for ( int i = 0; i < Val.cols(); i++ ) { |
308 | | x ( i ) = evallog ( Val.get_col ( i ) ) ; |
309 | | } |
310 | | return x; |
311 | | } |
312 | | //! Compute log-probability of multiple values argument \c val |
313 | | virtual vec evallog_m ( const Array<vec> &Avec ) const { |
314 | | vec x ( Avec.size() ); |
315 | | for ( int i = 0; i < Avec.size(); i++ ) { |
316 | | x ( i ) = evallog ( Avec ( i ) ) ; |
317 | | } |
318 | | return x; |
319 | | } |
320 | | //! Return conditional density on the given RV, the remaining rvs will be in conditioning |
321 | | virtual mpdf* condition ( const RV &rv ) const { |
322 | | it_warning ( "Not implemented" ); |
323 | | return NULL; |
324 | | } |
325 | | |
326 | | //! Return marginal density on the given RV, the remainig rvs are intergrated out |
327 | | virtual epdf* marginal ( const RV &rv ) const { |
328 | | it_warning ( "Not implemented" ); |
329 | | return NULL; |
330 | | } |
331 | | |
332 | | //! return expected value |
333 | | virtual vec mean() const { |
334 | | it_error ( "not implemneted" ); |
335 | | return vec ( 0 ); |
336 | | }; |
337 | | |
338 | | //! return expected variance (not covariance!) |
339 | | virtual vec variance() const { |
340 | | it_error ( "not implemneted" ); |
341 | | return vec ( 0 ); |
342 | | }; |
343 | | //! Lower and upper bounds of \c percentage % quantile, returns mean-2*sigma as default |
344 | | virtual void qbounds ( vec &lb, vec &ub, double percentage = 0.95 ) const { |
345 | | vec mea = mean(); |
346 | | vec std = sqrt ( variance() ); |
347 | | lb = mea - 2 * std; |
348 | | ub = mea + 2 * std; |
349 | | }; |
350 | | //!@} |
351 | | |
352 | | //! \name Connection to other classes |
353 | | //! Description of the random quantity via attribute \c rv is optional. |
354 | | //! For operations such as sampling \c rv does not need to be set. However, for \c marginalization |
355 | | //! and \c conditioning \c rv has to be set. NB: |
356 | | //! @{ |
357 | | |
358 | | //!Name its rv |
359 | | void set_rv ( const RV &rv0 ) { |
360 | | rv = rv0; |
361 | | } //it_assert_debug(isnamed(),""); }; |
362 | | //! True if rv is assigned |
363 | | bool isnamed() const { |
364 | | bool b = ( dim == rv._dsize() ); |
365 | | return b; |
366 | | } |
367 | | //! Return name (fails when isnamed is false) |
368 | | const RV& _rv() const { |
369 | | it_assert_debug ( isnamed(), "" ); |
370 | | return rv; |
371 | | } |
372 | | //!@} |
373 | | |
374 | | //! \name Access to attributes |
375 | | //! @{ |
376 | | |
377 | | //! Size of the random variable |
378 | | int dimension() const { |
379 | | return dim; |
380 | | } |
381 | | //! Load from structure with elements: |
382 | | //! \code |
383 | | //! { rv = {class="RV", names=(...),}; // RV describing meaning of random variable |
384 | | //! // elements of offsprings |
385 | | //! } |
386 | | //! \endcode |
387 | | //!@} |
388 | | void from_setting ( const Setting &set ) { |
389 | | RV* r = UI::build<RV> ( set, "rv" ); |
390 | | if ( r ) { |
391 | | set_rv ( *r ); |
392 | | delete r; |
393 | | } |
394 | | } |
| 263 | class epdf : public root |
| 264 | { |
| 265 | protected: |
| 266 | //! dimension of the random variable |
| 267 | int dim; |
| 268 | //! Description of the random variable |
| 269 | RV rv; |
| 270 | |
| 271 | public: |
| 272 | /*! \name Constructors |
| 273 | Construction of each epdf should support two types of constructors: |
| 274 | \li empty constructor, |
| 275 | \li copy constructor, |
| 276 | |
| 277 | The following constructors should be supported for convenience: |
| 278 | \li constructor followed by calling \c set_parameters() |
| 279 | \li constructor accepting random variables calling \c set_rv() |
| 280 | |
| 281 | All internal data structures are constructed as empty. Their values (including sizes) will be set by method \c set_parameters(). This way references can be initialized in constructors. |
| 282 | @{*/ |
| 283 | epdf() : dim (0), rv() {}; |
| 284 | epdf (const epdf &e) : dim (e.dim), rv (e.rv) {}; |
| 285 | epdf (const RV &rv0) : dim (rv0._dsize()) { |
| 286 | set_rv (rv0); |
| 287 | }; |
| 288 | void set_parameters (int dim0) { |
| 289 | dim = dim0; |
| 290 | } |
| 291 | //!@} |
| 292 | |
| 293 | //! \name Matematical Operations |
| 294 | //!@{ |
| 295 | |
| 296 | //! Returns a sample, \f$ x \f$ from density \f$ f_x()\f$ |
| 297 | virtual vec sample() const { |
| 298 | it_error ("not implemneted"); |
| 299 | return vec (0); |
| 300 | }; |
| 301 | //! Returns N samples, \f$ [x_1 , x_2 , \ldots \ \f$ from density \f$ f_x(rv)\f$ |
| 302 | virtual mat sample_m (int N) const; |
| 303 | //! Compute log-probability of argument \c val |
| 304 | //! In case the argument is out of suport return -Infinity |
| 305 | virtual double evallog (const vec &val) const { |
| 306 | it_error ("not implemneted"); |
| 307 | return 0.0; |
| 308 | }; |
| 309 | //! Compute log-probability of multiple values argument \c val |
| 310 | virtual vec evallog_m (const mat &Val) const { |
| 311 | vec x (Val.cols()); |
| 312 | for (int i = 0; i < Val.cols(); i++) { |
| 313 | x (i) = evallog (Val.get_col (i)) ; |
| 314 | } |
| 315 | return x; |
| 316 | } |
| 317 | //! Compute log-probability of multiple values argument \c val |
| 318 | virtual vec evallog_m (const Array<vec> &Avec) const { |
| 319 | vec x (Avec.size()); |
| 320 | for (int i = 0; i < Avec.size(); i++) { |
| 321 | x (i) = evallog (Avec (i)) ; |
| 322 | } |
| 323 | return x; |
| 324 | } |
| 325 | //! Return conditional density on the given RV, the remaining rvs will be in conditioning |
| 326 | virtual mpdf* condition (const RV &rv) const { |
| 327 | it_warning ("Not implemented"); |
| 328 | return NULL; |
| 329 | } |
| 330 | |
| 331 | //! Return marginal density on the given RV, the remainig rvs are intergrated out |
| 332 | virtual epdf* marginal (const RV &rv) const { |
| 333 | it_warning ("Not implemented"); |
| 334 | return NULL; |
| 335 | } |
| 336 | |
| 337 | //! return expected value |
| 338 | virtual vec mean() const { |
| 339 | it_error ("not implemneted"); |
| 340 | return vec (0); |
| 341 | }; |
| 342 | |
| 343 | //! return expected variance (not covariance!) |
| 344 | virtual vec variance() const { |
| 345 | it_error ("not implemneted"); |
| 346 | return vec (0); |
| 347 | }; |
| 348 | //! Lower and upper bounds of \c percentage % quantile, returns mean-2*sigma as default |
| 349 | virtual void qbounds (vec &lb, vec &ub, double percentage = 0.95) const { |
| 350 | vec mea = mean(); |
| 351 | vec std = sqrt (variance()); |
| 352 | lb = mea - 2 * std; |
| 353 | ub = mea + 2 * std; |
| 354 | }; |
| 355 | //!@} |
| 356 | |
| 357 | //! \name Connection to other classes |
| 358 | //! Description of the random quantity via attribute \c rv is optional. |
| 359 | //! For operations such as sampling \c rv does not need to be set. However, for \c marginalization |
| 360 | //! and \c conditioning \c rv has to be set. NB: |
| 361 | //! @{ |
| 362 | |
| 363 | //!Name its rv |
| 364 | void set_rv (const RV &rv0) { |
| 365 | rv = rv0; |
| 366 | } //it_assert_debug(isnamed(),""); }; |
| 367 | //! True if rv is assigned |
| 368 | bool isnamed() const { |
| 369 | bool b = (dim == rv._dsize()); |
| 370 | return b; |
| 371 | } |
| 372 | //! Return name (fails when isnamed is false) |
| 373 | const RV& _rv() const { |
| 374 | it_assert_debug (isnamed(), ""); |
| 375 | return rv; |
| 376 | } |
| 377 | //!@} |
| 378 | |
| 379 | //! \name Access to attributes |
| 380 | //! @{ |
| 381 | |
| 382 | //! Size of the random variable |
| 383 | int dimension() const { |
| 384 | return dim; |
| 385 | } |
| 386 | //! Load from structure with elements: |
| 387 | //! \code |
| 388 | //! { rv = {class="RV", names=(...),}; // RV describing meaning of random variable |
| 389 | //! // elements of offsprings |
| 390 | //! } |
| 391 | //! \endcode |
| 392 | //!@} |
| 393 | void from_setting (const Setting &set) { |
| 394 | RV* r = UI::build<RV> (set, "rv"); |
| 395 | if (r) { |
| 396 | set_rv (*r); |
| 397 | delete r; |
| 398 | } |
| 399 | } |
401 | | class mpdf : public root { |
402 | | protected: |
403 | | //!dimension of the condition |
404 | | int dimc; |
405 | | //! random variable in condition |
406 | | RV rvc; |
407 | | private: |
408 | | //! internal epdf, used only as cache to avoid virtual calls of \c _rv() and \c _dimension() |
409 | | epdf* ep; |
410 | | |
411 | | protected: |
412 | | void set_ep(epdf &iepdf) { |
413 | | ep = &iepdf; |
414 | | } |
415 | | |
416 | | public: |
417 | | //! \name Constructors |
418 | | //! @{ |
419 | | |
420 | | mpdf() : dimc ( 0 ), rvc(), ep(NULL) { } |
421 | | |
422 | | mpdf ( const mpdf &m ) : dimc ( m.dimc ), rvc ( m.rvc ), ep ( m.ep ) { } |
423 | | //!@} |
424 | | |
425 | | //! \name Matematical operations |
426 | | //!@{ |
427 | | |
428 | | //! Returns a sample from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv |
429 | | virtual vec samplecond ( const vec &cond ){it_error("Not implemented");return vec(0);}; |
430 | | |
431 | | //! Returns \param N samples from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv |
432 | | virtual mat samplecond_m ( const vec &cond, int N ){it_error("Not implemented");return mat();} |
433 | | |
434 | | |
435 | | //! Shortcut for conditioning and evaluation of the internal epdf. In some cases, this operation can be implemented efficiently. |
436 | | virtual double evallogcond ( const vec &dt, const vec &cond ){it_error("Not implemented");return 0.0;} |
437 | | |
438 | | //! Matrix version of evallogcond |
439 | | virtual vec evallogcond_m ( const mat &Dt, const vec &cond ){it_error("Not implemented");return vec(0);} |
440 | | |
441 | | //! Array<vec> version of evallogcond |
442 | | virtual vec evallogcond_m ( const Array<vec> &Dt, const vec &cond ){it_error("Not implemented");return vec(0);} |
443 | | |
444 | | //! \name Access to attributes |
445 | | //! @{ |
446 | | |
447 | | RV _rv() { |
448 | | return ep->_rv(); |
449 | | } |
450 | | RV _rvc() { |
451 | | it_assert_debug ( isnamed(), "" ); |
452 | | return rvc; |
453 | | } |
454 | | int dimension() { |
455 | | return ep->dimension(); |
456 | | } |
457 | | int dimensionc() { |
458 | | return dimc; |
459 | | } |
460 | | |
461 | | //! Load from structure with elements: |
462 | | //! \code |
463 | | //! { class = "mpdf_offspring", |
464 | | //! rv = {class="RV", names=(...),}; // RV describing meaning of random variable |
465 | | //! rvc= {class="RV", names=(...),}; // RV describing meaning of random variable in condition |
466 | | //! // elements of offsprings |
467 | | //! } |
468 | | //! \endcode |
469 | | //!@} |
470 | | void from_setting ( const Setting &set ); |
471 | | //!@} |
472 | | |
473 | | //! \name Connection to other objects |
474 | | //!@{ |
475 | | void set_rvc ( const RV &rvc0 ) { |
476 | | rvc = rvc0; |
477 | | } |
478 | | void set_rv ( const RV &rv0 ) { |
479 | | ep->set_rv ( rv0 ); |
480 | | } |
481 | | bool isnamed() { |
482 | | return ( ep->isnamed() ) && ( dimc == rvc._dsize() ); |
483 | | } |
484 | | //!@} |
| 406 | class mpdf : public root |
| 407 | { |
| 408 | protected: |
| 409 | //!dimension of the condition |
| 410 | int dimc; |
| 411 | //! random variable in condition |
| 412 | RV rvc; |
| 413 | private: |
| 414 | //! internal epdf, used only as cache to avoid virtual calls of \c _rv() and \c _dimension() |
| 415 | epdf* ep; |
| 416 | |
| 417 | protected: |
| 418 | void set_ep (epdf &iepdf) { |
| 419 | ep = &iepdf; |
| 420 | } |
| 421 | |
| 422 | public: |
| 423 | //! \name Constructors |
| 424 | //! @{ |
| 425 | |
| 426 | mpdf() : dimc (0), rvc(), ep (NULL) { } |
| 427 | |
| 428 | mpdf (const mpdf &m) : dimc (m.dimc), rvc (m.rvc), ep (m.ep) { } |
| 429 | //!@} |
| 430 | |
| 431 | //! \name Matematical operations |
| 432 | //!@{ |
| 433 | |
| 434 | //! Returns a sample from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv |
| 435 | virtual vec samplecond (const vec &cond) {it_error ("Not implemented");return vec (0);}; |
| 436 | |
| 437 | //! Returns \param N samples from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv |
| 438 | virtual mat samplecond_m (const vec &cond, int N) {it_error ("Not implemented");return mat();} |
| 439 | |
| 440 | |
| 441 | //! Shortcut for conditioning and evaluation of the internal epdf. In some cases, this operation can be implemented efficiently. |
| 442 | virtual double evallogcond (const vec &dt, const vec &cond) {it_error ("Not implemented");return 0.0;} |
| 443 | |
| 444 | //! Matrix version of evallogcond |
| 445 | virtual vec evallogcond_m (const mat &Dt, const vec &cond) {it_error ("Not implemented");return vec (0);} |
| 446 | |
| 447 | //! Array<vec> version of evallogcond |
| 448 | virtual vec evallogcond_m (const Array<vec> &Dt, const vec &cond) {it_error ("Not implemented");return vec (0);} |
| 449 | |
| 450 | //! \name Access to attributes |
| 451 | //! @{ |
| 452 | |
| 453 | RV _rv() { |
| 454 | return ep->_rv(); |
| 455 | } |
| 456 | RV _rvc() { |
| 457 | it_assert_debug (isnamed(), ""); |
| 458 | return rvc; |
| 459 | } |
| 460 | int dimension() { |
| 461 | return ep->dimension(); |
| 462 | } |
| 463 | int dimensionc() { |
| 464 | return dimc; |
| 465 | } |
| 466 | |
| 467 | //! Load from structure with elements: |
| 468 | //! \code |
| 469 | //! { class = "mpdf_offspring", |
| 470 | //! rv = {class="RV", names=(...),}; // RV describing meaning of random variable |
| 471 | //! rvc= {class="RV", names=(...),}; // RV describing meaning of random variable in condition |
| 472 | //! // elements of offsprings |
| 473 | //! } |
| 474 | //! \endcode |
| 475 | //!@} |
| 476 | void from_setting (const Setting &set); |
| 477 | //!@} |
| 478 | |
| 479 | //! \name Connection to other objects |
| 480 | //!@{ |
| 481 | void set_rvc (const RV &rvc0) { |
| 482 | rvc = rvc0; |
| 483 | } |
| 484 | void set_rv (const RV &rv0) { |
| 485 | ep->set_rv (rv0); |
| 486 | } |
| 487 | bool isnamed() { |
| 488 | return (ep->isnamed()) && (dimc == rvc._dsize()); |
| 489 | } |
| 490 | //!@} |
539 | | class datalink { |
540 | | protected: |
541 | | //! Remember how long val should be |
542 | | int downsize; |
543 | | |
544 | | //! Remember how long val of "Up" should be |
545 | | int upsize; |
546 | | |
547 | | //! val-to-val link, indices of the upper val |
548 | | ivec v2v_up; |
549 | | |
550 | | public: |
551 | | //! Constructor |
552 | | datalink() : downsize ( 0 ), upsize ( 0 ) { } |
553 | | datalink ( const RV &rv, const RV &rv_up ) { |
554 | | set_connection ( rv, rv_up ); |
555 | | } |
556 | | |
557 | | //! set connection, rv must be fully present in rv_up |
558 | | void set_connection ( const RV &rv, const RV &rv_up ) { |
559 | | downsize = rv._dsize(); |
560 | | upsize = rv_up._dsize(); |
561 | | v2v_up = rv.dataind ( rv_up ); |
562 | | |
563 | | it_assert_debug ( v2v_up.length() == downsize, "rv is not fully in rv_up" ); |
564 | | } |
565 | | |
566 | | //! set connection using indices |
567 | | void set_connection ( int ds, int us, const ivec &upind ) { |
568 | | downsize = ds; |
569 | | upsize = us; |
570 | | v2v_up = upind; |
571 | | |
572 | | it_assert_debug ( v2v_up.length() == downsize, "rv is not fully in rv_up" ); |
573 | | } |
574 | | |
575 | | //! Get val for myself from val of "Up" |
576 | | vec pushdown ( const vec &val_up ) { |
577 | | it_assert_debug ( upsize == val_up.length(), "Wrong val_up" ); |
578 | | return get_vec ( val_up, v2v_up ); |
579 | | } |
580 | | |
581 | | //! Fill val of "Up" by my pieces |
582 | | void pushup ( vec &val_up, const vec &val ) { |
583 | | it_assert_debug ( downsize == val.length(), "Wrong val" ); |
584 | | it_assert_debug ( upsize == val_up.length(), "Wrong val_up" ); |
585 | | set_subvector ( val_up, v2v_up, val ); |
586 | | } |
| 546 | class datalink |
| 547 | { |
| 548 | protected: |
| 549 | //! Remember how long val should be |
| 550 | int downsize; |
| 551 | |
| 552 | //! Remember how long val of "Up" should be |
| 553 | int upsize; |
| 554 | |
| 555 | //! val-to-val link, indices of the upper val |
| 556 | ivec v2v_up; |
| 557 | |
| 558 | public: |
| 559 | //! Constructor |
| 560 | datalink() : downsize (0), upsize (0) { } |
| 561 | datalink (const RV &rv, const RV &rv_up) { |
| 562 | set_connection (rv, rv_up); |
| 563 | } |
| 564 | |
| 565 | //! set connection, rv must be fully present in rv_up |
| 566 | void set_connection (const RV &rv, const RV &rv_up) { |
| 567 | downsize = rv._dsize(); |
| 568 | upsize = rv_up._dsize(); |
| 569 | v2v_up = rv.dataind (rv_up); |
| 570 | |
| 571 | it_assert_debug (v2v_up.length() == downsize, "rv is not fully in rv_up"); |
| 572 | } |
| 573 | |
| 574 | //! set connection using indices |
| 575 | void set_connection (int ds, int us, const ivec &upind) { |
| 576 | downsize = ds; |
| 577 | upsize = us; |
| 578 | v2v_up = upind; |
| 579 | |
| 580 | it_assert_debug (v2v_up.length() == downsize, "rv is not fully in rv_up"); |
| 581 | } |
| 582 | |
| 583 | //! Get val for myself from val of "Up" |
| 584 | vec pushdown (const vec &val_up) { |
| 585 | it_assert_debug (upsize == val_up.length(), "Wrong val_up"); |
| 586 | return get_vec (val_up, v2v_up); |
| 587 | } |
| 588 | |
| 589 | //! Fill val of "Up" by my pieces |
| 590 | void pushup (vec &val_up, const vec &val) { |
| 591 | it_assert_debug (downsize == val.length(), "Wrong val"); |
| 592 | it_assert_debug (upsize == val_up.length(), "Wrong val_up"); |
| 593 | set_subvector (val_up, v2v_up, val); |
| 594 | } |
590 | | class datalink_m2e: public datalink { |
591 | | protected: |
592 | | //! Remember how long cond should be |
593 | | int condsize; |
594 | | |
595 | | //!upper_val-to-local_cond link, indices of the upper val |
596 | | ivec v2c_up; |
597 | | |
598 | | //!upper_val-to-local_cond link, indices of the local cond |
599 | | ivec v2c_lo; |
600 | | |
601 | | public: |
602 | | //! Constructor |
603 | | datalink_m2e() : condsize ( 0 ) { } |
604 | | |
605 | | void set_connection ( const RV &rv, const RV &rvc, const RV &rv_up ) { |
606 | | datalink::set_connection ( rv, rv_up ); |
607 | | condsize = rvc._dsize(); |
608 | | //establish v2c connection |
609 | | rvc.dataind ( rv_up, v2c_lo, v2c_up ); |
610 | | } |
611 | | |
612 | | //!Construct condition |
613 | | vec get_cond ( const vec &val_up ) { |
614 | | vec tmp ( condsize ); |
615 | | set_subvector ( tmp, v2c_lo, val_up ( v2c_up ) ); |
616 | | return tmp; |
617 | | } |
618 | | |
619 | | void pushup_cond ( vec &val_up, const vec &val, const vec &cond ) { |
620 | | it_assert_debug ( downsize == val.length(), "Wrong val" ); |
621 | | it_assert_debug ( upsize == val_up.length(), "Wrong val_up" ); |
622 | | set_subvector ( val_up, v2v_up, val ); |
623 | | set_subvector ( val_up, v2c_up, cond ); |
624 | | } |
| 598 | class datalink_m2e: public datalink |
| 599 | { |
| 600 | protected: |
| 601 | //! Remember how long cond should be |
| 602 | int condsize; |
| 603 | |
| 604 | //!upper_val-to-local_cond link, indices of the upper val |
| 605 | ivec v2c_up; |
| 606 | |
| 607 | //!upper_val-to-local_cond link, indices of the local cond |
| 608 | ivec v2c_lo; |
| 609 | |
| 610 | public: |
| 611 | //! Constructor |
| 612 | datalink_m2e() : condsize (0) { } |
| 613 | |
| 614 | void set_connection (const RV &rv, const RV &rvc, const RV &rv_up) { |
| 615 | datalink::set_connection (rv, rv_up); |
| 616 | condsize = rvc._dsize(); |
| 617 | //establish v2c connection |
| 618 | rvc.dataind (rv_up, v2c_lo, v2c_up); |
| 619 | } |
| 620 | |
| 621 | //!Construct condition |
| 622 | vec get_cond (const vec &val_up) { |
| 623 | vec tmp (condsize); |
| 624 | set_subvector (tmp, v2c_lo, val_up (v2c_up)); |
| 625 | return tmp; |
| 626 | } |
| 627 | |
| 628 | void pushup_cond (vec &val_up, const vec &val, const vec &cond) { |
| 629 | it_assert_debug (downsize == val.length(), "Wrong val"); |
| 630 | it_assert_debug (upsize == val_up.length(), "Wrong val_up"); |
| 631 | set_subvector (val_up, v2v_up, val); |
| 632 | set_subvector (val_up, v2c_up, cond); |
| 633 | } |
629 | | class datalink_m2m: public datalink_m2e { |
630 | | protected: |
631 | | //!cond-to-cond link, indices of the upper cond |
632 | | ivec c2c_up; |
633 | | //!cond-to-cond link, indices of the local cond |
634 | | ivec c2c_lo; |
635 | | |
636 | | public: |
637 | | //! Constructor |
638 | | datalink_m2m() {}; |
639 | | void set_connection ( const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up ) { |
640 | | datalink_m2e::set_connection ( rv, rvc, rv_up ); |
641 | | //establish c2c connection |
642 | | rvc.dataind ( rvc_up, c2c_lo, c2c_up ); |
643 | | it_assert_debug ( c2c_lo.length() + v2c_lo.length() == condsize, "cond is not fully given" ); |
644 | | } |
645 | | |
646 | | //! Get cond for myself from val and cond of "Up" |
647 | | vec get_cond ( const vec &val_up, const vec &cond_up ) { |
648 | | vec tmp ( condsize ); |
649 | | set_subvector ( tmp, v2c_lo, val_up ( v2c_up ) ); |
650 | | set_subvector ( tmp, c2c_lo, cond_up ( c2c_up ) ); |
651 | | return tmp; |
652 | | } |
653 | | //! Fill |
| 638 | class datalink_m2m: public datalink_m2e |
| 639 | { |
| 640 | protected: |
| 641 | //!cond-to-cond link, indices of the upper cond |
| 642 | ivec c2c_up; |
| 643 | //!cond-to-cond link, indices of the local cond |
| 644 | ivec c2c_lo; |
| 645 | |
| 646 | public: |
| 647 | //! Constructor |
| 648 | datalink_m2m() {}; |
| 649 | void set_connection (const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up) { |
| 650 | datalink_m2e::set_connection (rv, rvc, rv_up); |
| 651 | //establish c2c connection |
| 652 | rvc.dataind (rvc_up, c2c_lo, c2c_up); |
| 653 | it_assert_debug (c2c_lo.length() + v2c_lo.length() == condsize, "cond is not fully given"); |
| 654 | } |
| 655 | |
| 656 | //! Get cond for myself from val and cond of "Up" |
| 657 | vec get_cond (const vec &val_up, const vec &cond_up) { |
| 658 | vec tmp (condsize); |
| 659 | set_subvector (tmp, v2c_lo, val_up (v2c_up)); |
| 660 | set_subvector (tmp, c2c_lo, cond_up (c2c_up)); |
| 661 | return tmp; |
| 662 | } |
| 663 | //! Fill |
662 | | class logger : public root { |
663 | | protected: |
664 | | //! RVs of all logged variables. |
665 | | Array<RV> entries; |
666 | | //! Names of logged quantities, e.g. names of algorithm variants |
667 | | Array<string> names; |
668 | | public: |
669 | | //!Default constructor |
670 | | logger() : entries ( 0 ), names ( 0 ) {} |
671 | | |
672 | | //! returns an identifier which will be later needed for calling the \c logit() function |
673 | | //! For empty RV it returns -1, this entry will be ignored by \c logit(). |
674 | | virtual int add ( const RV &rv, string prefix = "" ) { |
675 | | int id; |
676 | | if ( rv._dsize() > 0 ) { |
677 | | id = entries.length(); |
678 | | names = concat ( names, prefix ); // diff |
679 | | entries.set_length ( id + 1, true ); |
680 | | entries ( id ) = rv; |
681 | | } else { |
682 | | id = -1; |
683 | | } |
684 | | return id; // identifier of the last entry |
685 | | } |
686 | | |
687 | | //! log this vector |
688 | | virtual void logit ( int id, const vec &v ) = 0; |
689 | | //! log this double |
690 | | virtual void logit ( int id, const double &d ) = 0; |
691 | | |
692 | | //! Shifts storage position for another time step. |
693 | | virtual void step() = 0; |
694 | | |
695 | | //! Finalize storing information |
696 | | virtual void finalize() {}; |
697 | | |
698 | | //! Initialize the storage |
699 | | virtual void init() {}; |
| 672 | class logger : public root |
| 673 | { |
| 674 | protected: |
| 675 | //! RVs of all logged variables. |
| 676 | Array<RV> entries; |
| 677 | //! Names of logged quantities, e.g. names of algorithm variants |
| 678 | Array<string> names; |
| 679 | public: |
| 680 | //!Default constructor |
| 681 | logger() : entries (0), names (0) {} |
| 682 | |
| 683 | //! returns an identifier which will be later needed for calling the \c logit() function |
| 684 | //! For empty RV it returns -1, this entry will be ignored by \c logit(). |
| 685 | virtual int add (const RV &rv, string prefix = "") { |
| 686 | int id; |
| 687 | if (rv._dsize() > 0) { |
| 688 | id = entries.length(); |
| 689 | names = concat (names, prefix); // diff |
| 690 | entries.set_length (id + 1, true); |
| 691 | entries (id) = rv; |
| 692 | } else { |
| 693 | id = -1; |
| 694 | } |
| 695 | return id; // identifier of the last entry |
| 696 | } |
| 697 | |
| 698 | //! log this vector |
| 699 | virtual void logit (int id, const vec &v) = 0; |
| 700 | //! log this double |
| 701 | virtual void logit (int id, const double &d) = 0; |
| 702 | |
| 703 | //! Shifts storage position for another time step. |
| 704 | virtual void step() = 0; |
| 705 | |
| 706 | //! Finalize storing information |
| 707 | virtual void finalize() {}; |
| 708 | |
| 709 | //! Initialize the storage |
| 710 | virtual void init() {}; |
768 | | class DS : public root { |
769 | | protected: |
770 | | int dtsize; |
771 | | int utsize; |
772 | | //!Description of data returned by \c getdata(). |
773 | | RV Drv; |
774 | | //!Description of data witten by by \c write(). |
775 | | RV Urv; // |
776 | | //! Remember its own index in Logger L |
777 | | int L_dt, L_ut; |
778 | | public: |
779 | | //! default constructors |
780 | | DS() : Drv(), Urv() {}; |
781 | | //! Returns full vector of observed data=[output, input] |
782 | | virtual void getdata ( vec &dt ) { |
783 | | it_error ( "abstract class" ); |
784 | | }; |
785 | | //! Returns data records at indeces. |
786 | | virtual void getdata ( vec &dt, const ivec &indeces ) { |
787 | | it_error ( "abstract class" ); |
788 | | }; |
789 | | //! Accepts action variable and schedule it for application. |
790 | | virtual void write ( vec &ut ) { |
791 | | it_error ( "abstract class" ); |
792 | | }; |
793 | | //! Accepts action variables at specific indeces |
794 | | virtual void write ( vec &ut, const ivec &indeces ) { |
795 | | it_error ( "abstract class" ); |
796 | | }; |
797 | | |
798 | | //! Moves from \f$ t \f$ to \f$ t+1 \f$, i.e. perfroms the actions and reads response of the system. |
799 | | virtual void step() = 0; |
800 | | |
801 | | //! Register DS for logging into logger L |
802 | | virtual void log_add ( logger &L ) { |
803 | | it_assert_debug ( dtsize == Drv._dsize(), "" ); |
804 | | it_assert_debug ( utsize == Urv._dsize(), "" ); |
805 | | |
806 | | L_dt = L.add ( Drv, "" ); |
807 | | L_ut = L.add ( Urv, "" ); |
808 | | } |
809 | | //! Register DS for logging into logger L |
810 | | virtual void logit ( logger &L ) { |
811 | | vec tmp ( Drv._dsize() + Urv._dsize() ); |
812 | | getdata ( tmp ); |
813 | | // d is first in getdata |
814 | | L.logit ( L_dt, tmp.left ( Drv._dsize() ) ); |
815 | | // u follows after d in getdata |
816 | | L.logit ( L_ut, tmp.mid ( Drv._dsize(), Urv._dsize() ) ); |
817 | | } |
818 | | //!access function |
819 | | virtual RV _drv() const { |
820 | | return concat ( Drv, Urv ); |
821 | | } |
822 | | //!access function |
823 | | const RV& _urv() const { |
824 | | return Urv; |
825 | | } |
826 | | //! set random rvariables |
827 | | virtual void set_drv ( const RV &drv, const RV &urv ) { |
828 | | Drv = drv; |
829 | | Urv = urv; |
830 | | } |
| 781 | class DS : public root |
| 782 | { |
| 783 | protected: |
| 784 | int dtsize; |
| 785 | int utsize; |
| 786 | //!Description of data returned by \c getdata(). |
| 787 | RV Drv; |
| 788 | //!Description of data witten by by \c write(). |
| 789 | RV Urv; // |
| 790 | //! Remember its own index in Logger L |
| 791 | int L_dt, L_ut; |
| 792 | public: |
| 793 | //! default constructors |
| 794 | DS() : Drv(), Urv() {}; |
| 795 | //! Returns full vector of observed data=[output, input] |
| 796 | virtual void getdata (vec &dt) { |
| 797 | it_error ("abstract class"); |
| 798 | }; |
| 799 | //! Returns data records at indeces. |
| 800 | virtual void getdata (vec &dt, const ivec &indeces) { |
| 801 | it_error ("abstract class"); |
| 802 | }; |
| 803 | //! Accepts action variable and schedule it for application. |
| 804 | virtual void write (vec &ut) { |
| 805 | it_error ("abstract class"); |
| 806 | }; |
| 807 | //! Accepts action variables at specific indeces |
| 808 | virtual void write (vec &ut, const ivec &indeces) { |
| 809 | it_error ("abstract class"); |
| 810 | }; |
| 811 | |
| 812 | //! Moves from \f$ t \f$ to \f$ t+1 \f$, i.e. perfroms the actions and reads response of the system. |
| 813 | virtual void step() = 0; |
| 814 | |
| 815 | //! Register DS for logging into logger L |
| 816 | virtual void log_add (logger &L) { |
| 817 | it_assert_debug (dtsize == Drv._dsize(), ""); |
| 818 | it_assert_debug (utsize == Urv._dsize(), ""); |
| 819 | |
| 820 | L_dt = L.add (Drv, ""); |
| 821 | L_ut = L.add (Urv, ""); |
| 822 | } |
| 823 | //! Register DS for logging into logger L |
| 824 | virtual void logit (logger &L) { |
| 825 | vec tmp (Drv._dsize() + Urv._dsize()); |
| 826 | getdata (tmp); |
| 827 | // d is first in getdata |
| 828 | L.logit (L_dt, tmp.left (Drv._dsize())); |
| 829 | // u follows after d in getdata |
| 830 | L.logit (L_ut, tmp.mid (Drv._dsize(), Urv._dsize())); |
| 831 | } |
| 832 | //!access function |
| 833 | virtual RV _drv() const { |
| 834 | return concat (Drv, Urv); |
| 835 | } |
| 836 | //!access function |
| 837 | const RV& _urv() const { |
| 838 | return Urv; |
| 839 | } |
| 840 | //! set random rvariables |
| 841 | virtual void set_drv (const RV &drv, const RV &urv) { |
| 842 | Drv = drv; |
| 843 | Urv = urv; |
| 844 | } |
854 | | class BM : public root { |
855 | | protected: |
856 | | //! Random variable of the data (optional) |
857 | | RV drv; |
858 | | //!Logarithm of marginalized data likelihood. |
859 | | double ll; |
860 | | //! If true, the filter will compute likelihood of the data record and store it in \c ll . Set to false if you want to save computational time. |
861 | | bool evalll; |
862 | | public: |
863 | | //! \name Constructors |
864 | | //! @{ |
865 | | |
866 | | BM() : ll ( 0 ), evalll ( true ), LIDs ( 4 ), LFlags ( 4 ) { |
867 | | LIDs = -1;/*empty IDs*/ |
868 | | LFlags = 0; |
869 | | LFlags ( 0 ) = 1;/*log only mean*/ |
870 | | }; |
871 | | BM ( const BM &B ) : drv ( B.drv ), ll ( B.ll ), evalll ( B.evalll ) {} |
872 | | //! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually! |
873 | | //! Prototype: \code BM* _copy_() const {return new BM(*this);} \endcode |
874 | | virtual BM* _copy_() const { |
875 | | return NULL; |
876 | | }; |
877 | | //!@} |
878 | | |
879 | | //! \name Mathematical operations |
880 | | //!@{ |
881 | | |
882 | | /*! \brief Incremental Bayes rule |
883 | | @param dt vector of input data |
884 | | */ |
885 | | virtual void bayes ( const vec &dt ) = 0; |
886 | | //! Batch Bayes rule (columns of Dt are observations) |
887 | | virtual void bayesB ( const mat &Dt ); |
888 | | //! Evaluates predictive log-likelihood of the given data record |
889 | | //! I.e. marginal likelihood of the data with the posterior integrated out. |
890 | | virtual double logpred ( const vec &dt ) const { |
891 | | it_error ( "Not implemented" ); |
892 | | return 0.0; |
| 868 | class BM : public root |
| 869 | { |
| 870 | protected: |
| 871 | //! Random variable of the data (optional) |
| 872 | RV drv; |
| 873 | //!Logarithm of marginalized data likelihood. |
| 874 | double ll; |
| 875 | //! If true, the filter will compute likelihood of the data record and store it in \c ll . Set to false if you want to save computational time. |
| 876 | bool evalll; |
| 877 | public: |
| 878 | //! \name Constructors |
| 879 | //! @{ |
| 880 | |
| 881 | BM() : ll (0), evalll (true), LIDs (4), LFlags (4) { |
| 882 | LIDs = -1;/*empty IDs*/ |
| 883 | LFlags = 0; |
| 884 | LFlags (0) = 1; /*log only mean*/ |
| 885 | }; |
| 886 | BM (const BM &B) : drv (B.drv), ll (B.ll), evalll (B.evalll) {} |
| 887 | //! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually! |
| 888 | //! Prototype: \code BM* _copy_() const {return new BM(*this);} \endcode |
| 889 | virtual BM* _copy_() const { |
| 890 | return NULL; |
| 891 | }; |
| 892 | //!@} |
| 893 | |
| 894 | //! \name Mathematical operations |
| 895 | //!@{ |
| 896 | |
| 897 | /*! \brief Incremental Bayes rule |
| 898 | @param dt vector of input data |
| 899 | */ |
| 900 | virtual void bayes (const vec &dt) = 0; |
| 901 | //! Batch Bayes rule (columns of Dt are observations) |
| 902 | virtual void bayesB (const mat &Dt); |
| 903 | //! Evaluates predictive log-likelihood of the given data record |
| 904 | //! I.e. marginal likelihood of the data with the posterior integrated out. |
| 905 | virtual double logpred (const vec &dt) const { |
| 906 | it_error ("Not implemented"); |
| 907 | return 0.0; |
| 908 | } |
| 909 | //! Matrix version of logpred |
| 910 | vec logpred_m (const mat &dt) const { |
| 911 | vec tmp (dt.cols()); |
| 912 | for (int i = 0; i < dt.cols(); i++) { |
| 913 | tmp (i) = logpred (dt.get_col (i)); |
| 914 | } |
| 915 | return tmp; |
| 916 | } |
| 917 | |
| 918 | //!Constructs a predictive density \f$ f(d_{t+1} |d_{t}, \ldots d_{0}) \f$ |
| 919 | virtual epdf* epredictor() const { |
| 920 | it_error ("Not implemented"); |
| 921 | return NULL; |
| 922 | }; |
| 923 | //!Constructs a conditional density 1-step ahead predictor \f$ f(d_{t+1} |d_{t+h-1}, \ldots d_{t}) |
| 924 | virtual mpdf* predictor() const { |
| 925 | it_error ("Not implemented"); |
| 926 | return NULL; |
| 927 | }; |
| 928 | //!@} |
| 929 | |
| 930 | //! \name Extension to conditional BM |
| 931 | //! This extension is useful e.g. in Marginalized Particle Filter (\ref bdm::MPF). |
| 932 | //! Alternatively, it can be used for automated connection to DS when the condition is observed |
| 933 | //!@{ |
| 934 | |
| 935 | //! Name of extension variable |
| 936 | RV rvc; |
| 937 | //! access function |
| 938 | const RV& _rvc() const { |
| 939 | return rvc; |
| 940 | } |
| 941 | |
| 942 | //! Substitute \c val for \c rvc. |
| 943 | virtual void condition (const vec &val) { |
| 944 | it_error ("Not implemented!"); |
| 945 | }; |
| 946 | |
| 947 | //!@} |
| 948 | |
| 949 | |
| 950 | //! \name Access to attributes |
| 951 | //!@{ |
| 952 | |
| 953 | const RV& _drv() const { |
| 954 | return drv; |
| 955 | } |
| 956 | void set_drv (const RV &rv) { |
| 957 | drv = rv; |
| 958 | } |
| 959 | void set_rv (const RV &rv) { |
| 960 | const_cast<epdf&> (posterior()).set_rv (rv); |
| 961 | } |
| 962 | double _ll() const { |
| 963 | return ll; |
| 964 | } |
| 965 | void set_evalll (bool evl0) { |
| 966 | evalll = evl0; |
| 967 | } |
| 968 | virtual const epdf& posterior() const = 0; |
| 969 | virtual const epdf* _e() const = 0; |
| 970 | //!@} |
| 971 | |
| 972 | //! \name Logging of results |
| 973 | //!@{ |
| 974 | |
| 975 | //! Set boolean options from a string, recognized are: "logbounds,logll" |
| 976 | virtual void set_options (const string &opt) { |
| 977 | LFlags (0) = 1; |
| 978 | if (opt.find ("logbounds") != string::npos) { |
| 979 | LFlags (1) = 1; |
| 980 | LFlags (2) = 1; |
| 981 | } |
| 982 | if (opt.find ("logll") != string::npos) { |
| 983 | LFlags (3) = 1; |
| 984 | } |
| 985 | } |
| 986 | //! IDs of storages in loggers 4:[1=mean,2=lb,3=ub,4=ll] |
| 987 | ivec LIDs; |
| 988 | |
| 989 | //! Flags for logging - same size as LIDs, each entry correspond to the same in LIDs |
| 990 | ivec LFlags; |
| 991 | //! Add all logged variables to a logger |
| 992 | virtual void log_add (logger &L, const string &name = "") { |
| 993 | // internal |
| 994 | RV r; |
| 995 | if (posterior().isnamed()) { |
| 996 | r = posterior()._rv(); |
| 997 | } else { |
| 998 | r = RV ("est", posterior().dimension()); |
| 999 | }; |
| 1000 | |
| 1001 | // Add mean value |
| 1002 | if (LFlags (0)) LIDs (0) = L.add (r, name + "mean_"); |
| 1003 | if (LFlags (1)) LIDs (1) = L.add (r, name + "lb_"); |
| 1004 | if (LFlags (2)) LIDs (2) = L.add (r, name + "ub_"); |
| 1005 | if (LFlags (3)) LIDs (3) = L.add (RV ("ll", 1), name); //TODO: "local" RV |
| 1006 | } |
| 1007 | virtual void logit (logger &L) { |
| 1008 | L.logit (LIDs (0), posterior().mean()); |
| 1009 | if (LFlags (1) || LFlags (2)) { //if one of them is off, its LID==-1 and will not be stored |
| 1010 | vec ub, lb; |
| 1011 | posterior().qbounds (lb, ub); |
| 1012 | L.logit (LIDs (1), lb); |
| 1013 | L.logit (LIDs (2), ub); |
| 1014 | } |
| 1015 | if (LFlags (3)) L.logit (LIDs (3), ll); |
| 1016 | } |
| 1017 | //!@} |
| 1018 | }; |
| 1019 | |
| 1020 | template<class EPDF> |
| 1021 | vec mpdf_internal<EPDF>::samplecond (const vec &cond) |
| 1022 | { |
| 1023 | condition (cond); |
| 1024 | vec temp = iepdf.sample(); |
| 1025 | return temp; |
| 1026 | } |
| 1027 | |
| 1028 | template<class EPDF> |
| 1029 | mat mpdf_internal<EPDF>::samplecond_m (const vec &cond, int N) |
| 1030 | { |
| 1031 | condition (cond); |
| 1032 | mat temp (dimension(), N); |
| 1033 | vec smp (dimension()); |
| 1034 | for (int i = 0; i < N; i++) { |
| 1035 | smp = iepdf.sample(); |
| 1036 | temp.set_col (i, smp); |
894 | | //! Matrix version of logpred |
895 | | vec logpred_m ( const mat &dt ) const { |
896 | | vec tmp ( dt.cols() ); |
897 | | for ( int i = 0; i < dt.cols(); i++ ) { |
898 | | tmp ( i ) = logpred ( dt.get_col ( i ) ); |
899 | | } |
900 | | return tmp; |
901 | | } |
902 | | |
903 | | //!Constructs a predictive density \f$ f(d_{t+1} |d_{t}, \ldots d_{0}) \f$ |
904 | | virtual epdf* epredictor() const { |
905 | | it_error ( "Not implemented" ); |
906 | | return NULL; |
907 | | }; |
908 | | //!Constructs a conditional density 1-step ahead predictor \f$ f(d_{t+1} |d_{t+h-1}, \ldots d_{t}) |
909 | | virtual mpdf* predictor() const { |
910 | | it_error ( "Not implemented" ); |
911 | | return NULL; |
912 | | }; |
913 | | //!@} |
914 | | |
915 | | //! \name Extension to conditional BM |
916 | | //! This extension is useful e.g. in Marginalized Particle Filter (\ref bdm::MPF). |
917 | | //! Alternatively, it can be used for automated connection to DS when the condition is observed |
918 | | //!@{ |
919 | | |
920 | | //! Name of extension variable |
921 | | RV rvc; |
922 | | //! access function |
923 | | const RV& _rvc() const { |
924 | | return rvc; |
925 | | } |
926 | | |
927 | | //! Substitute \c val for \c rvc. |
928 | | virtual void condition ( const vec &val ) { |
929 | | it_error ( "Not implemented!" ); |
930 | | }; |
931 | | |
932 | | //!@} |
933 | | |
934 | | |
935 | | //! \name Access to attributes |
936 | | //!@{ |
937 | | |
938 | | const RV& _drv() const { |
939 | | return drv; |
940 | | } |
941 | | void set_drv ( const RV &rv ) { |
942 | | drv = rv; |
943 | | } |
944 | | void set_rv ( const RV &rv ) { |
945 | | const_cast<epdf&> ( posterior() ).set_rv ( rv ); |
946 | | } |
947 | | double _ll() const { |
948 | | return ll; |
949 | | } |
950 | | void set_evalll ( bool evl0 ) { |
951 | | evalll = evl0; |
952 | | } |
953 | | virtual const epdf& posterior() const = 0; |
954 | | virtual const epdf* _e() const = 0; |
955 | | //!@} |
956 | | |
957 | | //! \name Logging of results |
958 | | //!@{ |
959 | | |
960 | | //! Set boolean options from a string, recognized are: "logbounds,logll" |
961 | | virtual void set_options ( const string &opt ) { |
962 | | LFlags ( 0 ) = 1; |
963 | | if ( opt.find ( "logbounds" ) != string::npos ) { |
964 | | LFlags ( 1 ) = 1; |
965 | | LFlags ( 2 ) = 1; |
966 | | } |
967 | | if ( opt.find ( "logll" ) != string::npos ) { |
968 | | LFlags ( 3 ) = 1; |
969 | | } |
970 | | } |
971 | | //! IDs of storages in loggers 4:[1=mean,2=lb,3=ub,4=ll] |
972 | | ivec LIDs; |
973 | | |
974 | | //! Flags for logging - same size as LIDs, each entry correspond to the same in LIDs |
975 | | ivec LFlags; |
976 | | //! Add all logged variables to a logger |
977 | | virtual void log_add ( logger &L, const string &name = "" ) { |
978 | | // internal |
979 | | RV r; |
980 | | if ( posterior().isnamed() ) { |
981 | | r = posterior()._rv(); |
982 | | } else { |
983 | | r = RV ( "est", posterior().dimension() ); |
984 | | }; |
985 | | |
986 | | // Add mean value |
987 | | if ( LFlags ( 0 ) ) LIDs ( 0 ) = L.add ( r, name + "mean_" ); |
988 | | if ( LFlags ( 1 ) ) LIDs ( 1 ) = L.add ( r, name + "lb_" ); |
989 | | if ( LFlags ( 2 ) ) LIDs ( 2 ) = L.add ( r, name + "ub_" ); |
990 | | if ( LFlags ( 3 ) ) LIDs ( 3 ) = L.add ( RV ( "ll", 1 ), name ); //TODO: "local" RV |
991 | | } |
992 | | virtual void logit ( logger &L ) { |
993 | | L.logit ( LIDs ( 0 ), posterior().mean() ); |
994 | | if ( LFlags ( 1 ) || LFlags ( 2 ) ) { //if one of them is off, its LID==-1 and will not be stored |
995 | | vec ub, lb; |
996 | | posterior().qbounds ( lb, ub ); |
997 | | L.logit ( LIDs ( 1 ), lb ); |
998 | | L.logit ( LIDs ( 2 ), ub ); |
999 | | } |
1000 | | if ( LFlags ( 3 ) ) L.logit ( LIDs ( 3 ), ll ); |
1001 | | } |
1002 | | //!@} |
1003 | | }; |
1004 | | |
| 1038 | |
| 1039 | return temp; |
| 1040 | } |
| 1041 | |
| 1042 | template<class EPDF> |
| 1043 | double mpdf_internal<EPDF>::evallogcond (const vec &dt, const vec &cond) |
| 1044 | { |
| 1045 | double tmp; |
| 1046 | condition (cond); |
| 1047 | tmp = iepdf.evallog (dt); |
| 1048 | // it_assert_debug(std::isfinite(tmp), "Infinite value"); |
| 1049 | return tmp; |
| 1050 | } |
| 1051 | |
| 1052 | template<class EPDF> |
| 1053 | vec mpdf_internal<EPDF>::evallogcond_m (const mat &Dt, const vec &cond) |
| 1054 | { |
| 1055 | condition (cond); |
| 1056 | return iepdf.evallog_m (Dt); |
| 1057 | } |
| 1058 | |
| 1059 | template<class EPDF> |
| 1060 | vec mpdf_internal<EPDF>::evallogcond_m (const Array<vec> &Dt, const vec &cond) |
| 1061 | { |
| 1062 | condition (cond); |
| 1063 | return iepdf.evallog_m (Dt); |
| 1064 | } |