| 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 | }; |