root/bdm/stat/libBM.h @ 214

/*!
  \file
  \brief Bayesian Models (bm) that use Bayes rule to learn from observations
  \author Vaclav Smidl.

  -----------------------------------
  BDM++ - C++ library for Bayesian Decision Making under Uncertainty

  Using IT++ for numerical operations
  -----------------------------------
*/

#ifndef BM_H
#define BM_H

#include <itpp/itbase.h>
#include "../itpp_ext.h"
//#include <std>

using namespace itpp;

//! Structure of RV (used internally), i.e. expanded RVs
class str {
public:
        //! vector id ids (non-unique!)
        ivec ids;
        //! vector of times
        ivec times;
        //!Default constructor
        str ( ivec ids0, ivec times0 ) :ids ( ids0 ),times ( times0 ) {
                it_assert_debug ( times0.length() ==ids0.length(),"Incompatible input" );
        };
};

/*!
* \brief Class representing variables, most often random variables

* More?...
*/

class RV {
protected:
        //! size = sum of sizes
        int tsize;
        //! len = number of individual rvs
        int len;
        //! Vector of unique IDs
        ivec ids;
        //! Vector of sizes
        ivec sizes;
        //! Vector of shifts from current time
        ivec times;
        //! Array of names
        Array<std::string> names;

private:
        //! auxiliary function used in constructor
        void init ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times );
public:
        //! Full constructor
        RV ( Array<std::string> in_names, ivec in_sizes, ivec in_times );
        //! Constructor with times=0
        RV ( Array<std::string> in_names, ivec in_sizes );
        //! Constructor with sizes=1, times=0
        RV ( Array<std::string> in_names );
        //! Constructor of empty RV
        RV ();

        //! Printing output e.g. for debugging.
        friend std::ostream &operator<< ( std::ostream &os, const RV &rv );

        //! Return number of scalars in the RV.
        int count() const {return tsize;} ;
        //! Return length (number of entries) of the RV.
        int length() const {return len;} ;

        //TODO why not inline and later??

        //! Find indexes of self in another rv, \return ivec of the same size as self.
        ivec findself ( const RV &rv2 ) const;
        //! Compare if \c rv2 is identical to this \c RV
        bool equal ( const RV &rv2 ) const;
        //! Add (concat) another variable to the current one, \return true if all rv2 were added, false if rv2 is in conflict
        bool add ( const RV &rv2 );
        //! Subtract  another variable from the current one
        RV subt ( const RV &rv2 ) const;
        //! Select only variables at indeces ind
        RV subselect ( const ivec &ind ) const;
        //! Select only variables at indeces ind
        RV operator() ( const ivec &ind ) const;
        //! Shift \c time shifted by delta.
        void t ( int delta );
        //! generate \c str from rv, by expanding sizes
        str tostr() const;
        //! when this rv is a part of bigger rv, this function returns indeces of self in the data vector of the bigger crv.
        //! Then, data can be copied via: data_of_this = cdata(ind);
        ivec dataind ( const RV &crv ) const;
        //! generate mutual indeces when copying data betwenn self and crv.
        //! Data are copied via: data_of_this(selfi) = data_of_rv2(rv2i)
        void dataind ( const RV &rv2, ivec &selfi, ivec &rv2i ) const;

        //!access function
        Array<std::string>& _names() {return names;};

        //!access function
        int id ( int at ) {return ids ( at );};
        //!access function
        int size ( int at ) {return sizes ( at );};
        //!access function
        int time ( int at ) {return times ( at );};
        //!access function
        std::string name ( int at ) {return names ( at );};

        //!access function
        void set_id ( int at, int id0 ) {ids ( at ) =id0;};
        //!access function
        void set_size ( int at, int size0 ) {sizes ( at ) =size0; tsize=sum ( sizes );};
        //!access function
        void set_time ( int at, int time0 ) {times ( at ) =time0;};

        //!Assign unused ids to this rv
        void newids();
};

//! Concat two random variables
RV concat ( const RV &rv1, const RV &rv2 );

//!Default empty RV that can be used as default argument
extern RV RV0;

//! Class representing function \f$f(x)\f$ of variable \f$x\f$ represented by \c rv

class fnc {
protected:
        //! Length of the output vector
        int dimy;
public:
        //!default constructor
        fnc ( int dy ) :dimy ( dy ) {};
        //! function evaluates numerical value of \f$f(x)\f$ at \f$x=\f$ \c cond
        virtual vec eval ( const vec &cond ) {
                return vec ( 0 );
        };

        //! access function
        int _dimy() const{return dimy;}

        //! Destructor for future use;
        virtual ~fnc() {};
};

class mpdf;

//! Probability density function with numerical statistics, e.g. posterior density.

class epdf {
protected:
        //! Identified of the random variable
        RV rv;
public:
        //!default constructor
        epdf() :rv ( ) {};

        //!default constructor
        epdf ( const RV &rv0 ) :rv ( rv0 ) {};

//      //! Returns the required moment of the epdf
//      virtual vec moment ( const int order = 1 );

        //! Returns a sample, \f$x\f$ from density \f$epdf(rv)\f$
        virtual vec sample () const =0;
        //! Returns N samples from density \f$epdf(rv)\f$
        virtual mat sample_m ( int N ) const;
        
        //! Compute log-probability of argument \c val
        virtual double evallog ( const vec &val ) const =0;

        //! Compute log-probability of multiple values argument \c val
        virtual vec evallog_m ( const mat &Val ) const {
                vec x ( Val.cols() );
                for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog ( Val.get_col ( i ) ) ;}
                return x;
        }
        //! Return conditional density on the given RV, the remaining rvs will be in conditioning
        virtual mpdf* condition ( const RV &rv ) const  {it_warning ( "Not implemented" ); return NULL;}
        //! Return marginal density on the given RV, the remainig rvs are intergrated out
        virtual epdf* marginal ( const RV &rv ) const {it_warning ( "Not implemented" ); return NULL;}

        //! return expected value
        virtual vec mean() const =0;

        //! Destructor for future use;
        virtual ~epdf() {};
        //! access function, possibly dangerous!
        const RV& _rv() const {return rv;}
        //! modifier function - useful when copying epdfs
        void _renewrv ( const RV &in_rv ) {rv=in_rv;}
        //!
};


//! Conditional probability density, e.g. modeling some dependencies.
//TODO Samplecond can be generalized

class mpdf {
protected:
        //! modeled random variable
        RV rv;
        //! random variable in condition
        RV rvc;
        //! pointer to internal epdf
        epdf* ep;
public:

        //! Returns the required moment of the epdf
//      virtual fnc moment ( const int order = 1 );
        //! 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.
        virtual vec samplecond ( const vec &cond, double &ll ) {
                this->condition ( cond );
                vec temp= ep->sample();
                ll=ep->evallog ( temp );return temp;
        };
        //! 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.
        virtual mat samplecond_m ( const vec &cond, vec &ll, int N ) {
                this->condition ( cond );
                mat temp ( rv.count(),N ); vec smp ( rv.count() );
                for ( int i=0;i<N;i++ ) {smp=ep->sample() ;temp.set_col ( i, smp );ll ( i ) =ep->evallog ( smp );}
                return temp;
        };
        //! Update \c ep so that it represents this mpdf conditioned on \c rvc = cond
        virtual void condition ( const vec &cond ) {it_error ( "Not implemented" );};

        //! Shortcut for conditioning and evaluation of the internal epdf. In some cases,  this operation can be implemented efficiently.
        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;
        };

        //! Matrix version of evallogcond
        virtual vec evallogcond_m ( const mat &Dt, const vec &cond ) {this->condition ( cond );return ep->evallog_m ( Dt );};

        //! Destructor for future use;
        virtual ~mpdf() {};

        //! Default constructor
        mpdf ( const RV &rv0, const RV &rvc0 ) :rv ( rv0 ),rvc ( rvc0 ) {};
        //! access function
        RV _rvc() const {return rvc;}
        //! access function
        RV _rv() const {return rv;}
        //!access function
        epdf& _epdf() {return *ep;}
};

//!DataLink is a connection between an epdf and its superordinate epdf (Up)
//! It is assumed that my val is fully present in "Up"
class datalink_e2e {
protected:
        //! Remember how long val should be
        int valsize;
        //! Remember how long val of "Up" should be
        int valupsize;
        //! val-to-val link, indeces of the upper val
        ivec v2v_up;
public:
        //! Constructor
        datalink_e2e ( const RV &rv, const RV &rv_up ) :
                        valsize ( rv.count() ), valupsize ( rv_up.count() ), v2v_up ( rv.dataind ( rv_up ) )  {
                it_assert_debug ( v2v_up.length() ==valsize,"rv is not fully in rv_up" );
        }
        //! Get val for myself from val of "Up"
        vec get_val ( const vec &val_up ) {it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" ); return get_vec ( val_up,v2v_up );}
        //! Fill val of "Up" by my pieces
        void fill_val ( vec &val_up, const vec &val ) {
                it_assert_debug ( valsize==val.length(),"Wrong val" );
                it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" );
                set_subvector ( val_up, v2v_up, val );
        }
};

//! data link between
class datalink_m2e: public datalink_e2e {
protected:
        //! Remember how long cond should be
        int condsize;
        //!upper_val-to-local_cond link, indeces of the upper val
        ivec v2c_up;
        //!upper_val-to-local_cond link, ideces of the local cond
        ivec v2c_lo;

public:
        //! Constructor
        datalink_m2e ( const RV &rv,  const RV &rvc, const RV &rv_up ) :
                        datalink_e2e ( rv,rv_up ), condsize ( rvc.count() ) {
                //establish v2c connection
                rvc.dataind ( rv_up, v2c_lo, v2c_up );
        }
        //!Construct condition
        vec get_cond ( const vec &val_up ) {
                vec tmp ( condsize );
                set_subvector ( tmp,v2c_lo,val_up ( v2c_up ) );
                return tmp;
        }
        void fill_val_cond ( vec &val_up, const vec &val, const vec &cond ) {
                it_assert_debug ( valsize==val.length(),"Wrong val" );
                it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" );
                set_subvector ( val_up, v2v_up, val );
                set_subvector ( val_up, v2c_up, cond );
        }
};
//!DataLink is a connection between mpdf and its superordinate (Up)
//! This class links
class datalink_m2m: public datalink_m2e {
protected:
        //!cond-to-cond link, indeces of the upper cond
        ivec c2c_up;
        //!cond-to-cond link, indeces of the local cond
        ivec c2c_lo;
public:
        //! Constructor
        datalink_m2m ( const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up ) :
                        datalink_m2e ( rv, rvc, rv_up) {
                //establish c2c connection
                rvc.dataind ( rvc_up, c2c_lo, c2c_up );
                it_assert_debug(c2c_lo.length()+v2c_lo.length()==condsize, "cond is not fully given");
        }
        //! Get cond for myself from val and cond of "Up"
        vec get_cond ( const vec &val_up, const vec &cond_up ) {
                vec tmp ( condsize );
                set_subvector ( tmp,v2c_lo,val_up ( v2c_up ) );
                set_subvector ( tmp,c2c_lo,cond_up ( c2c_up ) );
                return tmp;
        }
        //! Fill

};

/*! \brief Unconditional mpdf, allows using epdf in the role of mpdf.

*/
class mepdf : public mpdf {
public:
        //!Default constructor
        mepdf (const epdf* em ) :mpdf ( em->_rv(),RV() ) {ep=const_cast<epdf*>(em);};
        void condition ( const vec &cond ) {}
};

//!\brief Abstract composition of pdfs, a base for specific classes
//!this abstract class is common to epdf and mpdf
class compositepdf {
protected:
        //!Number of mpdfs in the composite
        int n;
        //! Elements of composition
        Array<mpdf*> mpdfs;
public:
        compositepdf ( Array<mpdf*> A0 ) : n ( A0.length() ), mpdfs ( A0 ) {};
        //! find common rv, flag \param checkoverlap modifies whether overlaps are acceptable
        RV getrv ( bool checkoverlap=false );
        //! common rvc of all mpdfs is written to rvc
        void setrvc ( const RV &rv, RV &rvc );
};

/*! \brief Abstract class for discrete-time sources of data.

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.
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).

*/

class DS {
protected:
        //!Observed variables, returned by \c getdata().
        RV Drv;
        //!Action variables, accepted by \c write().
        RV Urv; //
public:
        //! Returns full vector of observed data
        void getdata ( vec &dt );
        //! Returns data records at indeces.
        void getdata ( vec &dt, ivec &indeces );
        //! Accepts action variable and schedule it for application.
        void write ( vec &ut );
        //! Accepts action variables at specific indeces
        void write ( vec &ut, ivec &indeces );
        /*! \brief Method that assigns random variables to the datasource.
        Typically, the datasource will be constructed without knowledge of random variables. This method will associate existing variables with RVs.

        (Inherited from m3k, may be deprecated soon).
        */
        void linkrvs ( RV &drv, RV &urv );

        //! Moves from \f$t\f$ to \f$t+1\f$, i.e. perfroms the actions and reads response of the system.
        void step();

};

/*! \brief Bayesian Model of the world, i.e. all uncertainty is modeled by probabilities.

*/

class BM {
protected:
        //!Random variable of the posterior
        RV rv;
        //!Logarithm of marginalized data likelihood.
        double ll;
        //!  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.
        bool evalll;
public:

        //!Default constructor
        BM ( const RV &rv0, double ll0=0,bool evalll0=true ) :rv ( rv0 ), ll ( ll0 ),evalll ( evalll0 ) {//Fixme: test rv
        };
        //!Copy constructor
        BM ( const BM &B ) : rv ( B.rv ), ll ( B.ll ), evalll ( B.evalll ) {}

        /*! \brief Incremental Bayes rule
        @param dt vector of input data
        */
        virtual void bayes ( const vec &dt ) = 0;
        //! Batch Bayes rule (columns of Dt are observations)
        virtual void bayesB ( const mat &Dt );
        //! Returns a reference to the epdf representing posterior density on parameters. 
        virtual const epdf& _epdf() const =0;

        //! Returns a pointer to the epdf representing posterior density on parameters. Use with care!
        virtual const epdf* _e() const =0;

        //! Evaluates predictive log-likelihood of the given data record
        //! I.e. marginal likelihood of the data with the posterior integrated out.
        virtual double logpred ( const vec &dt ) const{it_error ( "Not implemented" );return 0.0;}
        //! Matrix version of logpred
        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;}

        //!Constructs a predictive density (marginal density on data)
        virtual epdf* predictor ( const RV &rv ) const {it_error ( "Not implemented" );return NULL;};

        //! Destructor for future use;
        virtual ~BM() {};
        //!access function
        const RV& _rv() const {return rv;}
        //!access function
        double _ll() const {return ll;}
        //!access function
        void set_evalll ( bool evl0 ) {evalll=evl0;}

        //! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually!
        //! Prototype: BM* _copy_(){BM Tmp*=new Tmp(this*);  return Tmp; }
        virtual BM* _copy_ ( bool changerv=false ) {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
};

/*!
\brief Conditional Bayesian Filter

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.

This is an interface class used to assure that certain BM has operation \c condition .

*/

class BMcond {
protected:
        //! Identificator of the conditioning variable
        RV rvc;
public:
        //! Substitute \c val for \c rvc.
        virtual void condition ( const vec &val ) =0;
        //! Default constructor
        BMcond ( RV &rv0 ) :rvc ( rv0 ) {};
        //! Destructor for future use
        virtual ~BMcond() {};
        //! access function
        const RV& _rvc() const {return rvc;}
};

#endif // BM_H