33 | | class str { |
34 | | public: |
35 | | //! vector id ids (non-unique!) |
36 | | ivec ids; |
37 | | //! vector of times |
38 | | ivec times; |
39 | | //!Default constructor |
40 | | str ( ivec ids0, ivec times0 ) :ids ( ids0 ),times ( times0 ) { |
41 | | it_assert_debug ( times0.length() ==ids0.length(),"Incompatible input" ); |
42 | | }; |
43 | | }; |
44 | | |
45 | | /*! |
46 | | * \brief Class representing variables, most often random variables |
47 | | |
48 | | * More?... |
49 | | */ |
50 | | |
51 | | class RV :public bdmroot{ |
52 | | protected: |
53 | | //! size = sum of sizes |
54 | | int tsize; |
55 | | //! len = number of individual rvs |
56 | | int len; |
57 | | //! Vector of unique IDs |
58 | | ivec ids; |
59 | | //! Vector of sizes |
60 | | ivec sizes; |
61 | | //! Vector of shifts from current time |
62 | | ivec times; |
63 | | //! Array of names |
64 | | Array<std::string> names; |
65 | | |
66 | | private: |
67 | | //! auxiliary function used in constructor |
68 | | void init ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times ); |
69 | | public: |
70 | | //! Full constructor |
71 | | RV ( Array<std::string> in_names, ivec in_sizes, ivec in_times ); |
72 | | //! Constructor with times=0 |
73 | | RV ( Array<std::string> in_names, ivec in_sizes ); |
74 | | //! Constructor with sizes=1, times=0 |
75 | | RV ( Array<std::string> in_names ); |
76 | | //! Constructor of empty RV |
77 | | RV (); |
78 | | |
79 | | //! Printing output e.g. for debugging. |
80 | | friend std::ostream &operator<< ( std::ostream &os, const RV &rv ); |
81 | | |
82 | | //! Return number of scalars in the RV. |
83 | | int count() const {return tsize;} ; |
84 | | //! Return length (number of entries) of the RV. |
85 | | int length() const {return len;} ; |
86 | | |
87 | | //TODO why not inline and later?? |
88 | | |
89 | | //! Find indexes of self in another rv, \return ivec of the same size as self. |
90 | | ivec findself ( const RV &rv2 ) const; |
91 | | //! Compare if \c rv2 is identical to this \c RV |
92 | | bool equal ( const RV &rv2 ) const; |
93 | | //! Add (concat) another variable to the current one, \return true if all rv2 were added, false if rv2 is in conflict |
94 | | bool add ( const RV &rv2 ); |
95 | | //! Subtract another variable from the current one |
96 | | RV subt ( const RV &rv2 ) const; |
97 | | //! Select only variables at indeces ind |
98 | | RV subselect ( const ivec &ind ) const; |
99 | | //! Select only variables at indeces ind |
100 | | RV operator() ( const ivec &ind ) const; |
101 | | //! Shift \c time shifted by delta. |
102 | | void t ( int delta ); |
103 | | //! generate \c str from rv, by expanding sizes |
104 | | str tostr() const; |
105 | | //! when this rv is a part of bigger rv, this function returns indeces of self in the data vector of the bigger crv. |
106 | | //! Then, data can be copied via: data_of_this = cdata(ind); |
107 | | ivec dataind ( const RV &crv ) const; |
108 | | //! generate mutual indeces when copying data betwenn self and crv. |
109 | | //! Data are copied via: data_of_this(selfi) = data_of_rv2(rv2i) |
110 | | void dataind ( const RV &rv2, ivec &selfi, ivec &rv2i ) const; |
111 | | |
112 | | //!access function |
113 | | Array<std::string>& _names() {return names;}; |
114 | | |
115 | | //!access function |
116 | | int id ( int at ) {return ids ( at );}; |
117 | | //!access function |
118 | | int size ( int at ) {return sizes ( at );}; |
119 | | //!access function |
120 | | int time ( int at ) {return times ( at );}; |
121 | | //!access function |
122 | | std::string name ( int at ) {return names ( at );}; |
123 | | |
124 | | //!access function |
125 | | void set_id ( int at, int id0 ) {ids ( at ) =id0;}; |
126 | | //!access function |
127 | | void set_size ( int at, int size0 ) {sizes ( at ) =size0; tsize=sum ( sizes );}; |
128 | | //!access function |
129 | | void set_time ( int at, int time0 ) {times ( at ) =time0;}; |
130 | | |
131 | | //!Assign unused ids to this rv |
132 | | void newids(); |
133 | | }; |
| 34 | class str { |
| 35 | public: |
| 36 | //! vector id ids (non-unique!) |
| 37 | ivec ids; |
| 38 | //! vector of times |
| 39 | ivec times; |
| 40 | //!Default constructor |
| 41 | str ( ivec ids0, ivec times0 ) :ids ( ids0 ),times ( times0 ) { |
| 42 | it_assert_debug ( times0.length() ==ids0.length(),"Incompatible input" ); |
| 43 | }; |
| 44 | }; |
| 45 | |
| 46 | /*! |
| 47 | * \brief Class representing variables, most often random variables |
| 48 | |
| 49 | * More?... |
| 50 | */ |
| 51 | |
| 52 | class RV :public bdmroot { |
| 53 | protected: |
| 54 | //! size = sum of sizes |
| 55 | int tsize; |
| 56 | //! len = number of individual rvs |
| 57 | int len; |
| 58 | //! Vector of unique IDs |
| 59 | ivec ids; |
| 60 | //! Vector of sizes |
| 61 | ivec sizes; |
| 62 | //! Vector of shifts from current time |
| 63 | ivec times; |
| 64 | //! Array of names |
| 65 | Array<std::string> names; |
| 66 | |
| 67 | private: |
| 68 | //! auxiliary function used in constructor |
| 69 | void init ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times ); |
| 70 | public: |
| 71 | //! Full constructor |
| 72 | RV ( Array<std::string> in_names, ivec in_sizes, ivec in_times ); |
| 73 | //! Constructor with times=0 |
| 74 | RV ( Array<std::string> in_names, ivec in_sizes ); |
| 75 | //! Constructor with sizes=1, times=0 |
| 76 | RV ( Array<std::string> in_names ); |
| 77 | //! Constructor of empty RV |
| 78 | RV (); |
| 79 | //! Constructor of a single RV with given id |
| 80 | RV (string name, int id); |
| 81 | |
| 82 | //! Printing output e.g. for debugging. |
| 83 | friend std::ostream &operator<< ( std::ostream &os, const RV &rv ); |
| 84 | |
| 85 | //! Return number of scalars in the RV. |
| 86 | int count() const {return tsize;} ; |
| 87 | //! Return length (number of entries) of the RV. |
| 88 | int length() const {return len;} ; |
| 89 | |
| 90 | //TODO why not inline and later?? |
| 91 | |
| 92 | //! Find indexes of self in another rv, \return ivec of the same size as self. |
| 93 | ivec findself ( const RV &rv2 ) const; |
| 94 | //! Compare if \c rv2 is identical to this \c RV |
| 95 | bool equal ( const RV &rv2 ) const; |
| 96 | //! Add (concat) another variable to the current one, \return true if all rv2 were added, false if rv2 is in conflict |
| 97 | bool add ( const RV &rv2 ); |
| 98 | //! Subtract another variable from the current one |
| 99 | RV subt ( const RV &rv2 ) const; |
| 100 | //! Select only variables at indeces ind |
| 101 | RV subselect ( const ivec &ind ) const; |
| 102 | //! Select only variables at indeces ind |
| 103 | RV operator() ( const ivec &ind ) const; |
| 104 | //! Shift \c time shifted by delta. |
| 105 | void t ( int delta ); |
| 106 | //! generate \c str from rv, by expanding sizes |
| 107 | str tostr() const; |
| 108 | //! when this rv is a part of bigger rv, this function returns indeces of self in the data vector of the bigger crv. |
| 109 | //! Then, data can be copied via: data_of_this = cdata(ind); |
| 110 | ivec dataind ( const RV &crv ) const; |
| 111 | //! generate mutual indeces when copying data betwenn self and crv. |
| 112 | //! Data are copied via: data_of_this(selfi) = data_of_rv2(rv2i) |
| 113 | void dataind ( const RV &rv2, ivec &selfi, ivec &rv2i ) const; |
| 114 | //! Minimum time-offset |
| 115 | int mint () const {return min ( times );}; |
| 116 | |
| 117 | //!access function |
| 118 | Array<std::string>& _names() {return names;}; |
| 119 | |
| 120 | //!access function |
| 121 | int id ( int at ) {return ids ( at );}; |
| 122 | //!access function |
| 123 | int size ( int at ) {return sizes ( at );}; |
| 124 | //!access function |
| 125 | int time ( int at ) {return times ( at );}; |
| 126 | //!access function |
| 127 | std::string name ( int at ) {return names ( at );}; |
| 128 | |
| 129 | //!access function |
| 130 | void set_id ( int at, int id0 ) {ids ( at ) =id0;}; |
| 131 | //!access function |
| 132 | void set_size ( int at, int size0 ) {sizes ( at ) =size0; tsize=sum ( sizes );}; |
| 133 | //!access function |
| 134 | void set_time ( int at, int time0 ) {times ( at ) =time0;}; |
| 135 | |
| 136 | //!Assign unused ids to this rv |
| 137 | void newids(); |
| 138 | }; |
183 | | //! Returns a sample, \f$x\f$ from density \f$epdf(rv)\f$ |
184 | | virtual vec sample () const =0; |
185 | | //! Returns N samples from density \f$epdf(rv)\f$ |
186 | | virtual mat sample_m ( int N ) const; |
187 | | |
188 | | //! Compute log-probability of argument \c val |
189 | | virtual double evallog ( const vec &val ) const =0; |
190 | | |
191 | | //! Compute log-probability of multiple values argument \c val |
192 | | virtual vec evallog_m ( const mat &Val ) const { |
193 | | vec x ( Val.cols() ); |
194 | | for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog ( Val.get_col ( i ) ) ;} |
195 | | return x; |
196 | | } |
197 | | //! Return conditional density on the given RV, the remaining rvs will be in conditioning |
198 | | virtual mpdf* condition ( const RV &rv ) const {it_warning ( "Not implemented" ); return NULL;} |
199 | | //! Return marginal density on the given RV, the remainig rvs are intergrated out |
200 | | virtual epdf* marginal ( const RV &rv ) const {it_warning ( "Not implemented" ); return NULL;} |
201 | | |
202 | | //! return expected value |
203 | | virtual vec mean() const =0; |
204 | | |
205 | | //! return expected variance (not covariance!) |
206 | | virtual vec variance() const = 0; |
207 | | |
208 | | //! Destructor for future use; |
209 | | virtual ~epdf() {}; |
210 | | //! access function, possibly dangerous! |
211 | | const RV& _rv() const {return rv;} |
212 | | //! modifier function - useful when copying epdfs |
213 | | void _renewrv ( const RV &in_rv ) {rv=in_rv;} |
214 | | //! |
215 | | }; |
| 188 | //! Returns a sample, \f$x\f$ from density \f$epdf(rv)\f$ |
| 189 | virtual vec sample () const =0; |
| 190 | //! Returns N samples from density \f$epdf(rv)\f$ |
| 191 | virtual mat sample_m ( int N ) const; |
| 192 | |
| 193 | //! Compute log-probability of argument \c val |
| 194 | virtual double evallog ( const vec &val ) const =0; |
| 195 | |
| 196 | //! Compute log-probability of multiple values argument \c val |
| 197 | virtual vec evallog_m ( const mat &Val ) const { |
| 198 | vec x ( Val.cols() ); |
| 199 | for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog ( Val.get_col ( i ) ) ;} |
| 200 | return x; |
| 201 | } |
| 202 | //! Return conditional density on the given RV, the remaining rvs will be in conditioning |
| 203 | virtual mpdf* condition ( const RV &rv ) const {it_warning ( "Not implemented" ); return NULL;} |
| 204 | //! Return marginal density on the given RV, the remainig rvs are intergrated out |
| 205 | virtual epdf* marginal ( const RV &rv ) const {it_warning ( "Not implemented" ); return NULL;} |
| 206 | |
| 207 | //! return expected value |
| 208 | virtual vec mean() const =0; |
| 209 | |
| 210 | //! return expected variance (not covariance!) |
| 211 | virtual vec variance() const = 0; |
| 212 | |
| 213 | //! Destructor for future use; |
| 214 | virtual ~epdf() {}; |
| 215 | //! access function, possibly dangerous! |
| 216 | const RV& _rv() const {return rv;} |
| 217 | //! modifier function - useful when copying epdfs |
| 218 | void _renewrv ( const RV &in_rv ) {rv=in_rv;} |
| 219 | //! |
| 220 | }; |
221 | | class mpdf : public bdmroot{ |
222 | | protected: |
223 | | //! modeled random variable |
224 | | RV rv; |
225 | | //! random variable in condition |
226 | | RV rvc; |
227 | | //! pointer to internal epdf |
228 | | epdf* ep; |
229 | | public: |
230 | | |
231 | | //! 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 \param ll is a return value of log-likelihood of the sample. |
232 | | virtual vec samplecond ( const vec &cond, double &ll ) { |
233 | | this->condition ( cond ); |
234 | | vec temp= ep->sample(); |
235 | | ll=ep->evallog ( temp );return temp; |
236 | | }; |
237 | | //! 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 \param ll is a return value of log-likelihood of the sample. |
238 | | virtual mat samplecond_m ( const vec &cond, vec &ll, int N ) { |
239 | | this->condition ( cond ); |
240 | | mat temp ( rv.count(),N ); vec smp ( rv.count() ); |
241 | | for ( int i=0;i<N;i++ ) {smp=ep->sample() ;temp.set_col ( i, smp );ll ( i ) =ep->evallog ( smp );} |
242 | | return temp; |
243 | | }; |
244 | | //! Update \c ep so that it represents this mpdf conditioned on \c rvc = cond |
245 | | virtual void condition ( const vec &cond ) {it_error ( "Not implemented" );}; |
246 | | |
247 | | //! Shortcut for conditioning and evaluation of the internal epdf. In some cases, this operation can be implemented efficiently. |
248 | | virtual double evallogcond ( const vec &dt, const vec &cond ) {double tmp; this->condition ( cond );tmp = ep->evallog ( dt ); it_assert_debug(std::isfinite(tmp),"Infinite value"); return tmp; |
249 | | }; |
250 | | |
251 | | //! Matrix version of evallogcond |
252 | | virtual vec evallogcond_m ( const mat &Dt, const vec &cond ) {this->condition ( cond );return ep->evallog_m ( Dt );}; |
253 | | |
254 | | //! Destructor for future use; |
255 | | virtual ~mpdf() {}; |
256 | | |
257 | | //! Default constructor |
258 | | mpdf ( const RV &rv0, const RV &rvc0 ) :rv ( rv0 ),rvc ( rvc0 ) {}; |
259 | | //! access function |
260 | | RV _rvc() const {return rvc;} |
261 | | //! access function |
262 | | RV _rv() const {return rv;} |
263 | | //!access function |
264 | | epdf& _epdf() {return *ep;} |
265 | | //!access function |
266 | | epdf* _e() {return ep;} |
267 | | }; |
268 | | |
269 | | //!DataLink is a connection between an epdf and its superordinate epdf (Up) |
270 | | //! It is assumed that my val is fully present in "Up" |
271 | | class datalink_e2e { |
272 | | protected: |
273 | | //! Remember how long val should be |
274 | | int valsize; |
275 | | //! Remember how long val of "Up" should be |
276 | | int valupsize; |
277 | | //! val-to-val link, indeces of the upper val |
278 | | ivec v2v_up; |
279 | | public: |
280 | | //! Constructor |
281 | | datalink_e2e ( const RV &rv, const RV &rv_up ) : |
282 | | valsize ( rv.count() ), valupsize ( rv_up.count() ), v2v_up ( rv.dataind ( rv_up ) ) { |
283 | | it_assert_debug ( v2v_up.length() ==valsize,"rv is not fully in rv_up" ); |
| 226 | class mpdf : public bdmroot { |
| 227 | protected: |
| 228 | //! modeled random variable |
| 229 | RV rv; |
| 230 | //! random variable in condition |
| 231 | RV rvc; |
| 232 | //! pointer to internal epdf |
| 233 | epdf* ep; |
| 234 | public: |
| 235 | |
| 236 | //! 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 \param ll is a return value of log-likelihood of the sample. |
| 237 | virtual vec samplecond ( const vec &cond, double &ll ) { |
| 238 | this->condition ( cond ); |
| 239 | vec temp= ep->sample(); |
| 240 | ll=ep->evallog ( temp );return temp; |
| 241 | }; |
| 242 | //! 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 \param ll is a return value of log-likelihood of the sample. |
| 243 | virtual mat samplecond_m ( const vec &cond, vec &ll, int N ) { |
| 244 | this->condition ( cond ); |
| 245 | mat temp ( rv.count(),N ); vec smp ( rv.count() ); |
| 246 | for ( int i=0;i<N;i++ ) {smp=ep->sample() ;temp.set_col ( i, smp );ll ( i ) =ep->evallog ( smp );} |
| 247 | return temp; |
| 248 | }; |
| 249 | //! Update \c ep so that it represents this mpdf conditioned on \c rvc = cond |
| 250 | virtual void condition ( const vec &cond ) {it_error ( "Not implemented" );}; |
| 251 | |
| 252 | //! Shortcut for conditioning and evaluation of the internal epdf. In some cases, this operation can be implemented efficiently. |
| 253 | virtual double evallogcond ( const vec &dt, const vec &cond ) { |
| 254 | double tmp; this->condition ( cond );tmp = ep->evallog ( dt ); it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); return tmp; |
| 255 | }; |
| 256 | |
| 257 | //! Matrix version of evallogcond |
| 258 | virtual vec evallogcond_m ( const mat &Dt, const vec &cond ) {this->condition ( cond );return ep->evallog_m ( Dt );}; |
| 259 | |
| 260 | //! Destructor for future use; |
| 261 | virtual ~mpdf() {}; |
| 262 | |
| 263 | //! Default constructor |
| 264 | mpdf ( const RV &rv0, const RV &rvc0 ) :rv ( rv0 ),rvc ( rvc0 ) {}; |
| 265 | //! access function |
| 266 | RV _rvc() const {return rvc;} |
| 267 | //! access function |
| 268 | RV _rv() const {return rv;} |
| 269 | //!access function |
| 270 | epdf& _epdf() {return *ep;} |
| 271 | //!access function |
| 272 | epdf* _e() {return ep;} |
| 273 | }; |
| 274 | |
| 275 | /*! \brief DataLink is a connection between two data vectors Up and Down |
| 276 | |
| 277 | Up can be longer than Down. Down must be fully present in Up (TODO optional) |
| 278 | See chart: |
| 279 | \dot |
| 280 | digraph datalink { |
| 281 | node [shape=record]; |
| 282 | subgraph cluster0 { |
| 283 | label = "Up"; |
| 284 | up [label="<1>|<2>|<3>|<4>|<5>"]; |
| 285 | color = "white" |
293 | | }; |
| 293 | up:1 -> down:1; |
| 294 | up:3 -> down:2; |
| 295 | up:5 -> down:3; |
| 296 | } |
| 297 | \enddot |
| 298 | |
| 299 | */ |
| 300 | |
| 301 | /*! |
| 302 | @brief Class for storing results (and semi-results) of an experiment |
| 303 | |
| 304 | This class abstracts logging of results from implementation. This class replaces direct logging of results (e.g. to files or to global variables) by calling methods of a logger. Specializations of this abstract class for specific storage method are designed. |
| 305 | */ |
| 306 | class logger : public bdmroot { |
| 307 | protected: |
| 308 | //! RVs of all logged variables. |
| 309 | Array<RV> entries; |
| 310 | //! Names of logged quantities, e.g. names of algorithm variants |
| 311 | Array<string> names; |
| 312 | public: |
| 313 | //!Default constructor |
| 314 | logger ( ) : entries ( 0 ),names ( 0 ) {} |
| 315 | |
| 316 | //! returns an identifier which will be later needed for calling the log() function |
| 317 | virtual int add ( const RV &rv, string name="" ) { |
| 318 | int id=entries.length(); |
| 319 | names=concat ( names, name ); // diff |
| 320 | entries.set_length ( id+1,true ); |
| 321 | entries ( id ) = rv; |
| 322 | return id; // identifier of the last entry |
| 323 | } |
| 324 | |
| 325 | //! log this vector |
| 326 | virtual void logit ( int id, const vec &v ) =0; |
| 327 | |
| 328 | //! Shifts storage position for another time step. |
| 329 | virtual void step() =0; |
| 330 | |
| 331 | //! Finalize storing information |
| 332 | virtual void finalize() {}; |
| 333 | |
| 334 | //! Initialize the storage |
| 335 | virtual void init() {}; |
| 336 | |
| 337 | //! for future use |
| 338 | virtual ~logger() {}; |
| 339 | }; |
| 340 | |
| 341 | class datalink_e2e { |
| 342 | protected: |
| 343 | //! Remember how long val should be |
| 344 | int valsize; |
| 345 | //! Remember how long val of "Up" should be |
| 346 | int valupsize; |
| 347 | //! val-to-val link, indeces of the upper val |
| 348 | ivec v2v_up; |
| 349 | public: |
| 350 | //! Constructor |
| 351 | datalink_e2e ( const RV &rv, const RV &rv_up ) : |
| 352 | valsize ( rv.count() ), valupsize ( rv_up.count() ), v2v_up ( rv.dataind ( rv_up ) ) { |
| 353 | it_assert_debug ( v2v_up.length() ==valsize,"rv is not fully in rv_up" ); |
| 354 | } |
| 355 | //! Get val for myself from val of "Up" |
| 356 | vec get_val ( const vec &val_up ) {it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" ); return get_vec ( val_up,v2v_up );} |
| 357 | //! Fill val of "Up" by my pieces |
| 358 | void fill_val ( vec &val_up, const vec &val ) { |
| 359 | it_assert_debug ( valsize==val.length(),"Wrong val" ); |
| 360 | it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" ); |
| 361 | set_subvector ( val_up, v2v_up, val ); |
| 362 | } |
| 363 | }; |
296 | | class datalink_m2e: public datalink_e2e { |
297 | | protected: |
298 | | //! Remember how long cond should be |
299 | | int condsize; |
300 | | //!upper_val-to-local_cond link, indeces of the upper val |
301 | | ivec v2c_up; |
302 | | //!upper_val-to-local_cond link, ideces of the local cond |
303 | | ivec v2c_lo; |
304 | | |
305 | | public: |
306 | | //! Constructor |
307 | | datalink_m2e ( const RV &rv, const RV &rvc, const RV &rv_up ) : |
308 | | datalink_e2e ( rv,rv_up ), condsize ( rvc.count() ) { |
309 | | //establish v2c connection |
310 | | rvc.dataind ( rv_up, v2c_lo, v2c_up ); |
311 | | } |
312 | | //!Construct condition |
313 | | vec get_cond ( const vec &val_up ) { |
314 | | vec tmp ( condsize ); |
315 | | set_subvector ( tmp,v2c_lo,val_up ( v2c_up ) ); |
316 | | return tmp; |
317 | | } |
318 | | void fill_val_cond ( vec &val_up, const vec &val, const vec &cond ) { |
319 | | it_assert_debug ( valsize==val.length(),"Wrong val" ); |
320 | | it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" ); |
321 | | set_subvector ( val_up, v2v_up, val ); |
322 | | set_subvector ( val_up, v2c_up, cond ); |
323 | | } |
324 | | }; |
| 366 | class datalink_m2e: public datalink_e2e { |
| 367 | protected: |
| 368 | //! Remember how long cond should be |
| 369 | int condsize; |
| 370 | //!upper_val-to-local_cond link, indeces of the upper val |
| 371 | ivec v2c_up; |
| 372 | //!upper_val-to-local_cond link, ideces of the local cond |
| 373 | ivec v2c_lo; |
| 374 | |
| 375 | public: |
| 376 | //! Constructor |
| 377 | datalink_m2e ( const RV &rv, const RV &rvc, const RV &rv_up ) : |
| 378 | datalink_e2e ( rv,rv_up ), condsize ( rvc.count() ) { |
| 379 | //establish v2c connection |
| 380 | rvc.dataind ( rv_up, v2c_lo, v2c_up ); |
| 381 | } |
| 382 | //!Construct condition |
| 383 | vec get_cond ( const vec &val_up ) { |
| 384 | vec tmp ( condsize ); |
| 385 | set_subvector ( tmp,v2c_lo,val_up ( v2c_up ) ); |
| 386 | return tmp; |
| 387 | } |
| 388 | void fill_val_cond ( vec &val_up, const vec &val, const vec &cond ) { |
| 389 | it_assert_debug ( valsize==val.length(),"Wrong val" ); |
| 390 | it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" ); |
| 391 | set_subvector ( val_up, v2v_up, val ); |
| 392 | set_subvector ( val_up, v2c_up, cond ); |
| 393 | } |
| 394 | }; |
327 | | class datalink_m2m: public datalink_m2e { |
328 | | protected: |
329 | | //!cond-to-cond link, indeces of the upper cond |
330 | | ivec c2c_up; |
331 | | //!cond-to-cond link, indeces of the local cond |
332 | | ivec c2c_lo; |
333 | | public: |
334 | | //! Constructor |
335 | | datalink_m2m ( const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up ) : |
336 | | datalink_m2e ( rv, rvc, rv_up) { |
337 | | //establish c2c connection |
338 | | rvc.dataind ( rvc_up, c2c_lo, c2c_up ); |
339 | | it_assert_debug(c2c_lo.length()+v2c_lo.length()==condsize, "cond is not fully given"); |
340 | | } |
341 | | //! Get cond for myself from val and cond of "Up" |
342 | | vec get_cond ( const vec &val_up, const vec &cond_up ) { |
343 | | vec tmp ( condsize ); |
344 | | set_subvector ( tmp,v2c_lo,val_up ( v2c_up ) ); |
345 | | set_subvector ( tmp,c2c_lo,cond_up ( c2c_up ) ); |
346 | | return tmp; |
347 | | } |
348 | | //! Fill |
349 | | |
350 | | }; |
351 | | |
352 | | /*! \brief Unconditional mpdf, allows using epdf in the role of mpdf. |
353 | | |
354 | | */ |
355 | | class mepdf : public mpdf { |
356 | | public: |
357 | | //!Default constructor |
358 | | mepdf (const epdf* em ) :mpdf ( em->_rv(),RV() ) {ep=const_cast<epdf*>(em);}; |
359 | | void condition ( const vec &cond ) {} |
360 | | }; |
| 397 | class datalink_m2m: public datalink_m2e { |
| 398 | protected: |
| 399 | //!cond-to-cond link, indeces of the upper cond |
| 400 | ivec c2c_up; |
| 401 | //!cond-to-cond link, indeces of the local cond |
| 402 | ivec c2c_lo; |
| 403 | public: |
| 404 | //! Constructor |
| 405 | datalink_m2m ( const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up ) : |
| 406 | datalink_m2e ( rv, rvc, rv_up ) { |
| 407 | //establish c2c connection |
| 408 | rvc.dataind ( rvc_up, c2c_lo, c2c_up ); |
| 409 | it_assert_debug ( c2c_lo.length() +v2c_lo.length() ==condsize, "cond is not fully given" ); |
| 410 | } |
| 411 | //! Get cond for myself from val and cond of "Up" |
| 412 | vec get_cond ( const vec &val_up, const vec &cond_up ) { |
| 413 | vec tmp ( condsize ); |
| 414 | set_subvector ( tmp,v2c_lo,val_up ( v2c_up ) ); |
| 415 | set_subvector ( tmp,c2c_lo,cond_up ( c2c_up ) ); |
| 416 | return tmp; |
| 417 | } |
| 418 | //! Fill |
| 419 | |
| 420 | }; |
| 421 | |
| 422 | /*! \brief Unconditional mpdf, allows using epdf in the role of mpdf. |
| 423 | |
| 424 | */ |
| 425 | class mepdf : public mpdf { |
| 426 | public: |
| 427 | //!Default constructor |
| 428 | mepdf ( const epdf* em ) :mpdf ( em->_rv(),RV() ) {ep=const_cast<epdf*> ( em );}; |
| 429 | void condition ( const vec &cond ) {} |
| 430 | }; |
364 | | class compositepdf { |
365 | | protected: |
366 | | //!Number of mpdfs in the composite |
367 | | int n; |
368 | | //! Elements of composition |
369 | | Array<mpdf*> mpdfs; |
370 | | public: |
371 | | compositepdf ( Array<mpdf*> A0 ) : n ( A0.length() ), mpdfs ( A0 ) {}; |
372 | | //! find common rv, flag \param checkoverlap modifies whether overlaps are acceptable |
373 | | RV getrv ( bool checkoverlap=false ); |
374 | | //! common rvc of all mpdfs is written to rvc |
375 | | void setrvc ( const RV &rv, RV &rvc ); |
376 | | }; |
377 | | |
378 | | /*! \brief Abstract class for discrete-time sources of data. |
379 | | |
380 | | The class abstracts operations of: (i) data aquisition, (ii) data-preprocessing, (iii) scaling of data, and (iv) data resampling from the task of estimation and control. |
381 | | Moreover, for controlled systems, it is able to receive the desired control action and perform it in the next step. (Or as soon as possible). |
382 | | |
383 | | */ |
384 | | |
385 | | class DS : public bdmroot{ |
386 | | protected: |
387 | | //!Observed variables, returned by \c getdata(). |
388 | | RV Drv; |
389 | | //!Action variables, accepted by \c write(). |
390 | | RV Urv; // |
391 | | public: |
392 | | DS():Drv(RV0),Urv(RV0) {}; |
393 | | //! Returns full vector of observed data |
394 | | void getdata ( vec &dt ); |
395 | | //! Returns data records at indeces. |
396 | | void getdata ( vec &dt, ivec &indeces ); |
397 | | //! Accepts action variable and schedule it for application. |
398 | | void write ( vec &ut ); |
399 | | //! Accepts action variables at specific indeces |
400 | | void write ( vec &ut, ivec &indeces ); |
401 | | /*! \brief Method that assigns random variables to the datasource. |
402 | | Typically, the datasource will be constructed without knowledge of random variables. This method will associate existing variables with RVs. |
403 | | |
404 | | (Inherited from m3k, may be deprecated soon). |
405 | | */ |
406 | | void linkrvs ( RV &drv, RV &urv ); |
407 | | |
408 | | //! Moves from \f$t\f$ to \f$t+1\f$, i.e. perfroms the actions and reads response of the system. |
409 | | void step(); |
410 | | |
411 | | }; |
412 | | |
413 | | /*! \brief Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities. |
414 | | |
415 | | */ |
416 | | |
417 | | class BM :public bdmroot{ |
418 | | protected: |
419 | | //!Random variable of the posterior |
420 | | RV rv; |
421 | | //!Logarithm of marginalized data likelihood. |
422 | | double ll; |
423 | | //! 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. |
424 | | bool evalll; |
425 | | public: |
426 | | |
427 | | //!Default constructor |
428 | | BM ( const RV &rv0, double ll0=0,bool evalll0=true ) :rv ( rv0 ), ll ( ll0 ),evalll ( evalll0 ) {//Fixme: test rv |
429 | | }; |
430 | | //!Copy constructor |
431 | | BM ( const BM &B ) : rv ( B.rv ), ll ( B.ll ), evalll ( B.evalll ) {} |
432 | | |
433 | | /*! \brief Incremental Bayes rule |
434 | | @param dt vector of input data |
435 | | */ |
436 | | virtual void bayes ( const vec &dt ) = 0; |
437 | | //! Batch Bayes rule (columns of Dt are observations) |
438 | | virtual void bayesB ( const mat &Dt ); |
439 | | //! Returns a reference to the epdf representing posterior density on parameters. |
440 | | virtual const epdf& _epdf() const =0; |
441 | | |
442 | | //! Returns a pointer to the epdf representing posterior density on parameters. Use with care! |
443 | | virtual const epdf* _e() const =0; |
444 | | |
445 | | //! Evaluates predictive log-likelihood of the given data record |
446 | | //! I.e. marginal likelihood of the data with the posterior integrated out. |
447 | | virtual double logpred ( const vec &dt ) const{it_error ( "Not implemented" );return 0.0;} |
448 | | //! Matrix version of logpred |
449 | | vec logpred_m ( const mat &dt ) const{vec tmp ( dt.cols() );for ( int i=0;i<dt.cols();i++ ) {tmp ( i ) =logpred ( dt.get_col ( i ) );}return tmp;} |
450 | | |
451 | | //!Constructs a predictive density (marginal density on data) |
452 | | virtual epdf* predictor ( const RV &rv ) const {it_error ( "Not implemented" );return NULL;}; |
453 | | |
454 | | //! Destructor for future use; |
455 | | virtual ~BM() {}; |
456 | | //!access function |
457 | | const RV& _rv() const {return rv;} |
458 | | //!access function |
459 | | double _ll() const {return ll;} |
460 | | //!access function |
461 | | void set_evalll ( bool evl0 ) {evalll=evl0;} |
462 | | |
463 | | //! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually! |
464 | | //! Prototype: BM* _copy_(){BM Tmp*=new Tmp(this*); return Tmp; } |
465 | | virtual BM* _copy_ ( bool changerv=false ) {it_error ( "function _copy_ not implemented for this BM" ); return NULL;}; |
466 | | }; |
467 | | |
468 | | /*! |
469 | | \brief Conditional Bayesian Filter |
470 | | |
471 | | Evaluates conditional filtering density \f$f(rv|rvc,data)\f$ for a given \c rvc which is specified in each step by calling function \c condition. |
472 | | |
473 | | This is an interface class used to assure that certain BM has operation \c condition . |
474 | | |
475 | | */ |
476 | | |
477 | | class BMcond :public bdmroot{ |
478 | | protected: |
479 | | //! Identificator of the conditioning variable |
480 | | RV rvc; |
481 | | public: |
482 | | //! Substitute \c val for \c rvc. |
483 | | virtual void condition ( const vec &val ) =0; |
484 | | //! Default constructor |
485 | | BMcond ( RV &rv0 ) :rvc ( rv0 ) {}; |
486 | | //! Destructor for future use |
487 | | virtual ~BMcond() {}; |
488 | | //! access function |
489 | | const RV& _rvc() const {return rvc;} |
490 | | }; |
| 434 | class compositepdf { |
| 435 | protected: |
| 436 | //!Number of mpdfs in the composite |
| 437 | int n; |
| 438 | //! Elements of composition |
| 439 | Array<mpdf*> mpdfs; |
| 440 | public: |
| 441 | compositepdf ( Array<mpdf*> A0 ) : n ( A0.length() ), mpdfs ( A0 ) {}; |
| 442 | //! find common rv, flag \param checkoverlap modifies whether overlaps are acceptable |
| 443 | RV getrv ( bool checkoverlap=false ); |
| 444 | //! common rvc of all mpdfs is written to rvc |
| 445 | void setrvc ( const RV &rv, RV &rvc ); |
| 446 | }; |
| 447 | |
| 448 | /*! \brief Abstract class for discrete-time sources of data. |
| 449 | |
| 450 | The class abstracts operations of: (i) data aquisition, (ii) data-preprocessing, (iii) scaling of data, and (iv) data resampling from the task of estimation and control. |
| 451 | Moreover, for controlled systems, it is able to receive the desired control action and perform it in the next step. (Or as soon as possible). |
| 452 | |
| 453 | */ |
| 454 | |
| 455 | class DS : public bdmroot { |
| 456 | protected: |
| 457 | //!Observed variables, returned by \c getdata(). |
| 458 | RV Drv; |
| 459 | //!Action variables, accepted by \c write(). |
| 460 | RV Urv; // |
| 461 | //! Remember its own index in Logger L |
| 462 | int L_dt, L_ut; |
| 463 | public: |
| 464 | DS() :Drv ( RV0 ),Urv ( RV0 ) {}; |
| 465 | DS ( const RV &Drv0, const RV &Urv0 ) :Drv ( Drv0 ),Urv ( Urv0 ) {}; |
| 466 | //! Returns full vector of observed data=[output, input] |
| 467 | virtual void getdata ( vec &dt ) {it_error ( "abstract class" );}; |
| 468 | //! Returns data records at indeces. |
| 469 | virtual void getdata ( vec &dt, const ivec &indeces ) {it_error ( "abstract class" );}; |
| 470 | //! Accepts action variable and schedule it for application. |
| 471 | virtual void write ( vec &ut ) {it_error ( "abstract class" );}; |
| 472 | //! Accepts action variables at specific indeces |
| 473 | virtual void write ( vec &ut, const ivec &indeces ) {it_error ( "abstract class" );}; |
| 474 | |
| 475 | //! Moves from \f$t\f$ to \f$t+1\f$, i.e. perfroms the actions and reads response of the system. |
| 476 | virtual void step() =0; |
| 477 | |
| 478 | //! Register DS for logging into logger L |
| 479 | virtual void log_add ( logger &L ) { |
| 480 | L_dt=L.add ( Drv,"" ); |
| 481 | L_ut=L.add ( Urv,"" ); |
| 482 | } |
| 483 | //! Register DS for logging into logger L |
| 484 | virtual void logit ( logger &L ) { |
| 485 | vec tmp(Drv.count()+Urv.count()); |
| 486 | getdata(tmp); |
| 487 | // d is first in getdata |
| 488 | L.logit ( L_dt,tmp.left ( Drv.count() ) ); |
| 489 | // u follows after d in getdata |
| 490 | L.logit ( L_ut,tmp.mid ( Drv.count(), Urv.count() ) ); |
| 491 | } |
| 492 | }; |
| 493 | |
| 494 | /*! \brief Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities. |
| 495 | |
| 496 | */ |
| 497 | |
| 498 | class BM :public bdmroot { |
| 499 | protected: |
| 500 | //!Random variable of the posterior |
| 501 | RV rv; |
| 502 | //!Logarithm of marginalized data likelihood. |
| 503 | double ll; |
| 504 | //! 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. |
| 505 | bool evalll; |
| 506 | public: |
| 507 | |
| 508 | //!Default constructor |
| 509 | BM ( const RV &rv0, double ll0=0,bool evalll0=true ) :rv ( rv0 ), ll ( ll0 ),evalll ( evalll0 ) {//Fixme: test rv |
| 510 | }; |
| 511 | //!Copy constructor |
| 512 | BM ( const BM &B ) : rv ( B.rv ), ll ( B.ll ), evalll ( B.evalll ) {} |
| 513 | |
| 514 | /*! \brief Incremental Bayes rule |
| 515 | @param dt vector of input data |
| 516 | */ |
| 517 | virtual void bayes ( const vec &dt ) = 0; |
| 518 | //! Batch Bayes rule (columns of Dt are observations) |
| 519 | virtual void bayesB ( const mat &Dt ); |
| 520 | //! Returns a reference to the epdf representing posterior density on parameters. |
| 521 | virtual const epdf& _epdf() const =0; |
| 522 | |
| 523 | //! Returns a pointer to the epdf representing posterior density on parameters. Use with care! |
| 524 | virtual const epdf* _e() const =0; |
| 525 | |
| 526 | //! Evaluates predictive log-likelihood of the given data record |
| 527 | //! I.e. marginal likelihood of the data with the posterior integrated out. |
| 528 | virtual double logpred ( const vec &dt ) const{it_error ( "Not implemented" );return 0.0;} |
| 529 | //! Matrix version of logpred |
| 530 | vec logpred_m ( const mat &dt ) const{vec tmp ( dt.cols() );for ( int i=0;i<dt.cols();i++ ) {tmp ( i ) =logpred ( dt.get_col ( i ) );}return tmp;} |
| 531 | |
| 532 | //!Constructs a predictive density (marginal density on data) |
| 533 | virtual epdf* predictor ( const RV &rv ) const {it_error ( "Not implemented" );return NULL;}; |
| 534 | |
| 535 | //! Destructor for future use; |
| 536 | virtual ~BM() {}; |
| 537 | //!access function |
| 538 | const RV& _rv() const {return rv;} |
| 539 | //!access function |
| 540 | double _ll() const {return ll;} |
| 541 | //!access function |
| 542 | void set_evalll ( bool evl0 ) {evalll=evl0;} |
| 543 | |
| 544 | //! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually! |
| 545 | //! Prototype: BM* _copy_(){BM Tmp*=new Tmp(this*); return Tmp; } |
| 546 | virtual BM* _copy_ ( bool changerv=false ) {it_error ( "function _copy_ not implemented for this BM" ); return NULL;}; |
| 547 | }; |
| 548 | |
| 549 | /*! |
| 550 | \brief Conditional Bayesian Filter |
| 551 | |
| 552 | Evaluates conditional filtering density \f$f(rv|rvc,data)\f$ for a given \c rvc which is specified in each step by calling function \c condition. |
| 553 | |
| 554 | This is an interface class used to assure that certain BM has operation \c condition . |
| 555 | |
| 556 | */ |
| 557 | |
| 558 | class BMcond :public bdmroot { |
| 559 | protected: |
| 560 | //! Identificator of the conditioning variable |
| 561 | RV rvc; |
| 562 | public: |
| 563 | //! Substitute \c val for \c rvc. |
| 564 | virtual void condition ( const vec &val ) =0; |
| 565 | //! Default constructor |
| 566 | BMcond ( RV &rv0 ) :rvc ( rv0 ) {}; |
| 567 | //! Destructor for future use |
| 568 | virtual ~BMcond() {}; |
| 569 | //! access function |
| 570 | const RV& _rvc() const {return rvc;} |
| 571 | }; |