root/library/bdm/stat/emix.h @ 1064

/*!
  \file
  \brief Probability distributions for Mixtures of pdfs
  \author Vaclav Smidl.

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

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

#ifndef EMIX_H
#define EMIX_H

#define LOG2  0.69314718055995

#include "../shared_ptr.h"
#include "exp_family.h"

namespace bdm {

//this comes first because it is used inside emix!

/*! \brief Class representing ratio of two densities
which arise e.g. by applying the Bayes rule.
It represents density in the form:
\f[
f(rv|rvc) = \frac{f(rv,rvc)}{f(rvc)}
\f]
where \f$ f(rvc) = \int f(rv,rvc) d\ rv \f$.

In particular this type of arise by conditioning of a mixture model.

At present the only supported operation is evallogcond().
 */
class mratio: public pdf {
protected:
    //! Nominator in the form of pdf
    const epdf* nom;

    //!Denominator in the form of epdf
    shared_ptr<epdf> den;

    //!flag for destructor
    bool destroynom;
    //!datalink between conditional and nom
    datalink_m2e dl;
public:
    //!Default constructor. By default, the given epdf is not copied!
    //! It is assumed that this function will be used only temporarily.
    mratio ( const epdf* nom0, const RV &rv, bool copy = false ) : pdf ( ), dl ( ) {
        // adjust rv and rvc

        set_rv ( rv );
        dim = rv._dsize();

        rvc = nom0->_rv().subt ( rv );
        dimc = rvc._dsize();

        //prepare data structures
        if ( copy ) {
            bdm_error ( "todo" );
            // destroynom = true;
        } else {
            nom = nom0;
            destroynom = false;
        }
        bdm_assert_debug ( rvc.length() > 0, "Makes no sense to use this object!" );

        // build denominator
        den = nom->marginal ( rvc );
        dl.set_connection ( rv, rvc, nom0->_rv() );
    };

    double evallogcond ( const vec &val, const vec &cond ) {
        double tmp;
        vec nom_val ( dimension() + dimc );
        dl.pushup_cond ( nom_val, val, cond );
        tmp = exp ( nom->evallog ( nom_val ) - den->evallog ( cond ) );
        return tmp;
    }

    //! 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
    virtual vec samplecond ( const vec &cond ) NOT_IMPLEMENTED(0);

    //! Object takes ownership of nom and will destroy it
    void ownnom() {
        destroynom = true;
    }
    //! Default destructor
    ~mratio() {
        if ( destroynom ) {
            delete nom;
        }
    }


private:
    // not implemented
    mratio ( const mratio & );
    mratio &operator= ( const mratio & );
};

class emix; //forward

/*! Base class (Interface) for mixtures
*/
class emix_base : public epdf {
protected:
    //! reference to vector of weights
    vec &w;
    //! function returning ith component
    virtual const epdf * component(const int &i) const=0;

    virtual int no_coms() const=0;

public:

    emix_base(vec &w0): w(w0) {}

    void validate ();

    vec sample() const;

    vec mean() const;

    vec variance() const;

    double evallog ( const vec &val ) const;

    vec evallog_mat ( const mat &Val ) const;

    //! Auxiliary function that returns pdflog for each component
    mat evallog_coms ( const mat &Val ) const;

    shared_ptr<epdf> marginal ( const RV &rv ) const;
    //! Update already existing marginal density  \c target
    void marginal ( const RV &rv, emix &target ) const;
    shared_ptr<pdf> condition ( const RV &rv ) const;

    //Access methods
    //! returns a reference to the internal weights. Use with Care!
    vec& _w() {
        return w;
    }

    const vec& _w() const {
        return w;
    }
    //!access
    const epdf* _com(int i) const {
        return component(i);
    }

};

/*!
* \brief Mixture of epdfs

Density function:
\f[
f(x) = \sum_{i=1}^{n} w_{i} f_i(x), \quad \sum_{i=1}^n w_i = 1.
\f]
where \f$f_i(x)\f$ is any density on random variable \f$x\f$, called \a component,

*/
class emix : public emix_base {
protected:
    //! weights of the components
    vec weights;

    //! Component (epdfs)
    Array<shared_ptr<epdf> > Coms;

public:
    //! Default constructor
    emix ( ) : emix_base ( weights) { }

    const epdf* component(const int &i) const {
        return Coms(i).get();
    }
    void validate();


    int no_coms() const {
        return Coms.length();
    }

    //! Load from structure with elements:
    //!  \code
    //! { class='emix';
    //!   pdfs = (..., ...);     // list of pdfs in the mixture
    //!   weights = ( 0.5, 0.5 ); // weights of pdfs in the mixture
    //! }
    //! \endcode
    //!@}
    void from_setting ( const Setting &set );
    void to_setting  (Setting  &set) const;

    void set_rv ( const RV &rv ) {
        epdf::set_rv ( rv );
        for ( int i = 0; i < no_coms(); i++ ) {
            Coms( i )->set_rv ( rv );
        }
    }

    Array<shared_ptr<epdf> >& _Coms ( ) {
        return Coms;
    }
};
SHAREDPTR ( emix );
UIREGISTER ( emix );


/*! \brief Chain rule decomposition of epdf

Probability density in the form of Chain-rule decomposition:
\[
f(x_1,x_2,x_3) = f(x_1|x_2,x_3)f(x_2,x_3)f(x_3)
\]
Note that
*/
class mprod: public pdf {
private:
    Array<shared_ptr<pdf> > pdfs;

    //! Data link for each pdfs
    Array<shared_ptr<datalink_m2m> > dls;

public:
    //! \brief Default constructor
    mprod() { }

    /*!\brief Constructor from list of mFacs
    */
    mprod ( const Array<shared_ptr<pdf> > &mFacs ) {
        set_elements ( mFacs );
    }
    //! Set internal \c pdfs from given values
    void set_elements ( const Array<shared_ptr<pdf> > &mFacs );

    double evallogcond ( const vec &val, const vec &cond );

    vec evallogcond_mat ( const mat &Dt, const vec &cond );

    vec evallogcond_mat ( const Array<vec> &Dt, const vec &cond );

    //TODO smarter...
    vec samplecond ( const vec &cond ) {
        //! Ugly hack to help to discover if mpfs are not in proper order. Correct solution = check that explicitely.
        vec smp = std::numeric_limits<double>::infinity() * ones ( dimension() );
        vec smpi;
        // Hard assumption here!!! We are going backwards, to assure that samples that are needed from smp are already generated!
        for ( int i = ( pdfs.length() - 1 ); i >= 0; i-- ) {
            // generate contribution of this pdf
            smpi = pdfs ( i )->samplecond ( dls ( i )->get_cond ( smp , cond ) );
            // copy contribution of this pdf into smp
            dls ( i )->pushup ( smp, smpi );
        }
        return smp;
    }

    //! Load from structure with elements:
    //!  \code
    //! { class='mprod';
    //!   pdfs = (..., ...);     // list of pdfs in the order of chain rule
    //! }
    //! \endcode
    //!@}
    void from_setting ( const Setting &set ) ;
    void        to_setting  (Setting  &set) const;


};
UIREGISTER ( mprod );
SHAREDPTR ( mprod );


//! Product of independent epdfs. For dependent pdfs, use mprod.
class eprod_base: public epdf {
protected:
    //! Array of indices
    Array<datalink*> dls;
    //! interface for a factor
public:
    virtual const epdf* factor(int i) const NOT_IMPLEMENTED(NULL);
    //!number of factors
    virtual const int no_factors() const NOT_IMPLEMENTED(0);
    //! Default constructor
    eprod_base () :  dls (0) {};
    //! Set internal
    vec mean() const;

    vec variance() const;

    vec sample() const;

    double evallog ( const vec &val ) const;

    //!Destructor
    ~eprod_base() {
        for ( int i = 0; i < dls.length(); i++ ) {
            delete dls ( i );
        }
    }
    void validate() {
        epdf::validate();
        dls.set_length ( no_factors() );

        bool independent = true;
        dim = 0;
        for ( int i = 0; i < no_factors(); i++ ) {
            independent = rv.add ( factor ( i )->_rv() );
            dim += factor ( i )->dimension();
            bdm_assert_debug ( independent, "eprod:: given components are not independent." );
        };

        //
        int cumdim = 0;
        int dimi = 0;
        int i;
        for ( i = 0; i < no_factors(); i++ ) {
            dls ( i ) = new datalink;
            if ( isnamed() ) { // rvs are complete
                dls ( i )->set_connection ( factor ( i )->_rv() , rv );
            } else { //rvs are not reliable
                dimi = factor ( i )->dimension();
                dls ( i )->set_connection ( dimi, dim, linspace ( cumdim, cumdim + dimi - 1 ) );
                cumdim += dimi;
            }
        }

    }
};

class eprod: public eprod_base {
protected:
    Array<shared_ptr<epdf> > factors;
public:
    const epdf* factor(int i) const {
        return factors(i).get();
    }
    const int no_factors() const {
        return factors.length();
    }
    void set_parameters ( const Array<shared_ptr<epdf> > &epdfs0) {
        factors = epdfs0;
    }
    void from_setting(const Setting &set) {
        UI::get(factors,set,"pdfs",UI::compulsory);
    }
};
UIREGISTER(eprod);

//! similar to eprod but used only internally -- factors are external pointers
class eprod_internal: public eprod_base {
protected:
    Array<epdf* > factors;
    const epdf* factor(int i) const {
        return factors(i);
    }
    const int no_factors() const {
        return factors.length();
    }
public:
    void set_parameters ( const Array<epdf *> &epdfs0) {
        factors = epdfs0;
    }
};

/*! \brief Mixture of pdfs with constant weights, all pdfs are of equal RV and RVC

*/
class mmix : public pdf {
protected:
    //! Component (pdfs)
    Array<shared_ptr<pdf> > Coms;
    //!weights of the components
    vec w;
public:
    //!Default constructor
    mmix() : Coms ( 0 ) { };

    double evallogcond ( const vec &dt, const vec &cond ) {
        double ll = 0.0;
        for ( int i = 0; i < Coms.length(); i++ ) {
            ll += Coms ( i )->evallogcond ( dt, cond );
        }
        return ll;
    }

    vec samplecond ( const vec &cond );

    //! Load from structure with elements:
    //!  \code
    //! { class='mmix';
    //!   pdfs = (..., ...);     // list of pdfs in the mixture
    //!   weights = ( 0.5, 0.5 ); // weights of pdfs in the mixture
    //! }
    //! \endcode
    //!@}
    void from_setting ( const Setting &set );
    void        to_setting  (Setting  &set) const;
    virtual void validate();
};
SHAREDPTR ( mmix );
UIREGISTER ( mmix );


//! Base class for all BM running as parallel update of internal BMs

class ProdBMBase : public BM {
protected :
    Array<vec_from_vec> bm_yt;
    Array<vec_from_2vec> bm_cond;
    class eprod_bm : public eprod_base {
        ProdBMBase & pb;
    public :
        eprod_bm(ProdBMBase &pb0): pb(pb0) {}
        const epdf* factor(int i ) const {
            return &(pb.bm(i)->posterior());
        }
        const int no_factors() const {
            return pb.no_bms();
        }
    } est;
public:
    ProdBMBase():est(*this) {}
    virtual BM* bm(int i) NOT_IMPLEMENTED(NULL);
    virtual int no_bms() const {
        return 0;
    }
    const epdf& posterior() const {
        return est;
    }
    void set_prior(const epdf *pri) {
        const eprod_base* ep=dynamic_cast<const eprod_base*>(pri);
        if (ep) {
            bdm_assert(ep->no_factors()!=no_bms() , "Given prior has "+ num2str(ep->no_factors()) + " while this ProdBM has "+
                       num2str(no_bms()) + "BMs");
            for (int i=0; i<no_bms(); i++) {
                bm(i)->set_prior(ep->factor(i));
            }
        }
    }

    void validate() {
        est.validate();
        BM::validate();
        // set links
        bm_yt.set_length(no_bms());
        bm_cond.set_length(no_bms());

        //

        for (int i=0; i<no_bms(); i++) {
            yrv.add(bm(i)->_yrv());
            rvc.add(bm(i)->_rvc());
        }
        rvc=rvc.subt(yrv);

        dimy = yrv._dsize();
        dimc = rvc._dsize();

        for (int i=0; i<no_bms(); i++) {
            bm_yt(i).connect(bm(i)->_yrv(), yrv);
            bm_cond(i).connect(bm(i)->_rvc(), yrv, rvc);
        }
    }
    void bayes(const vec &dt, const vec &cond) {
        ll=0;
        for(int i=0; i<no_bms(); i++) {
            bm_yt(i).update(dt);
            bm_cond(i).update(dt,cond);
            bm(i)->bayes(bm_yt(i), bm_cond(i));
        }
    }

};

class ProdBM: public ProdBMBase {
protected:
    Array<shared_ptr<BM> > BMs;
public:
    virtual BM* bm(int i) {
        return BMs(i).get();
    }
    virtual int no_bms() const {
        return BMs.length();
    }
    void from_setting(const Setting &set) {
        BM::from_setting(set);
        UI::get(BMs,set,"BMs");
    }
    void to_setting(Setting &set) const {
        BM::to_setting(set);
        UI::save(BMs,set,"BMs");
    }

};
UIREGISTER(ProdBM);

//! \brief class for on-line model selection
class ModelComparator: public ProdBM {
protected:
    multiBM weights;
public:
    void bayes(const vec &yt, const vec &cond) {
        vec w_nn(no_bms());
        for (int i=0; i<no_bms(); i++) {
            bm(i)->bayes(yt,cond);
            w_nn(i) += bm(i)->_ll();
        }
        vec w=exp(w_nn-max(w_nn));
        weights.bayes(w/sum(w));
    }
    void validate() {
        ProdBM::validate();
        weights.validate();
    }
    void from_setting(const Setting& set) {
        ProdBM::from_setting(set);
        UI::get(weights.frg, set, "frg",UI::optional);
    }
    void to_setting(Setting& set) const {
        ProdBM::to_setting(set);
        UI::save(weights.frg, set, "frg");
    }
};

}
#endif //MX_H