emix.h

Go to the documentation of this file.

#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!

class mratio: public mpdf {
protected:
        const epdf* nom;
        epdf* den;
        bool destroynom;
        datalink_m2e dl;
public:
        mratio ( const epdf* nom0, const RV &rv, bool copy=false ) :mpdf ( ), dl ( ) {
                // adjust rv and rvc
                rvc = nom0->_rv().subt ( rv );
                dimc = rvc._dsize();
                set_ep(shared_ptr<epdf>(new epdf));
                e()->set_parameters(rv._dsize());
                e()->set_rv(rv);
                
                //prepare data structures
                if ( copy ) {it_error ( "todo" ); destroynom=true; }
                else { nom = nom0; destroynom = false; }
                it_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 ( e()->dimension() + dimc );
                dl.pushup_cond ( nom_val,val,cond );
                tmp = exp ( nom->evallog ( nom_val ) - den->evallog ( cond ) );
                it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
                return tmp;
        }
        void ownnom() {destroynom=true;}
        ~mratio() {delete den; if ( destroynom ) {delete nom;}}
};

class emix : public epdf {
protected:
        vec w;
        Array<epdf*> Coms;
        bool destroyComs;
public:
        emix ( ) : epdf ( ) {};
        void set_parameters ( const vec &w, const Array<epdf*> &Coms, bool copy=false );

        vec sample() const;
        vec mean() const {
                int i; vec mu = zeros ( dim );
                for ( i = 0;i < w.length();i++ ) {mu += w ( i ) * Coms ( i )->mean(); }
                return mu;
        }
        vec variance() const {
                //non-central moment
                vec mom2 = zeros ( dim );
                for ( int i = 0;i < w.length();i++ ) {mom2 += w ( i ) * (Coms(i)->variance() + pow ( Coms ( i )->mean(),2 )); }
                //central moment
                return mom2-pow ( mean(),2 );
        }
        double evallog ( const vec &val ) const {
                int i;
                double sum = 0.0;
                for ( i = 0;i < w.length();i++ ) {sum += w ( i ) * exp ( Coms ( i )->evallog ( val ) );}
                if ( sum==0.0 ) {sum=std::numeric_limits<double>::epsilon();}
                double tmp=log ( sum );
                it_assert_debug ( std::isfinite ( tmp ),"Infinite" );
                return tmp;
        };
        vec evallog_m ( const mat &Val ) const {
                vec x=zeros ( Val.cols() );
                for ( int i = 0; i < w.length(); i++ ) {
                        x+= w ( i ) *exp ( Coms ( i )->evallog_m ( Val ) );
                }
                return log ( x );
        };
        mat evallog_M ( const mat &Val ) const {
                mat X ( w.length(), Val.cols() );
                for ( int i = 0; i < w.length(); i++ ) {
                        X.set_row ( i, w ( i ) *exp ( Coms ( i )->evallog_m ( Val ) ) );
                }
                return X;
        };

        emix* marginal ( const RV &rv ) const;
        mratio* condition ( const RV &rv ) const; //why not mratio!!

//Access methods
        vec& _w() {return w;}
        virtual ~emix() {if ( destroyComs ) {for ( int i=0;i<Coms.length();i++ ) {delete Coms ( i );}}}
        void ownComs() {destroyComs=true;}

        epdf* _Coms ( int i ) {return Coms ( i );}
        void set_rv(const RV &rv){
                epdf::set_rv(rv);
                for(int i=0;i<Coms.length();i++){Coms(i)->set_rv(rv);}
        }
};


class egiwmix : public egiw {
protected:
        vec w;
        Array<egiw*> Coms;
        bool destroyComs;
public:
        egiwmix ( ) : egiw ( ) {};

        void set_parameters ( const vec &w, const Array<egiw*> &Coms, bool copy=false );

        vec mean() const;

        vec sample() const;

        vec variance() const;

        // TODO!!! Defined to follow ANSI and/or for future development
        void mean_mat ( mat &M, mat&R ) const {};
        double evallog_nn ( const vec &val ) const {return 0;};
        double lognc () const {return 0;};
        emix* marginal ( const RV &rv ) const;

//Access methods
        vec& _w() {return w;}
        virtual ~egiwmix() {if ( destroyComs ) {for ( int i=0;i<Coms.length();i++ ) {delete Coms ( i );}}}
        void ownComs() {destroyComs=true;}

        egiw* _Coms ( int i ) {return Coms ( i );}

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

        egiw* approx();
};

class mprod: public compositepdf, public mpdf {
protected:
        Array<epdf*> epdfs;
        Array<datalink_m2m*> dls;

        shared_ptr<epdf> dummy;

public:
        mprod():dummy(new epdf()) { }
        mprod (Array<mpdf*> mFacs):
            dummy(new epdf()) {
            set_elements(mFacs);
        }

        void set_elements(Array<mpdf*> mFacs , bool own=false) {
                
                compositepdf::set_elements(mFacs,own);
                dls.set_size(mFacs.length());
                epdfs.set_size(mFacs.length());
                                
                set_ep(dummy);
                RV rv=getrv ( true );
                set_rv ( rv );
                dummy->set_parameters(rv._dsize());
                setrvc(e()->_rv(), rvc);
                // rv and rvc established = > we can link them with mpdfs
                for ( int i = 0;i < mpdfs.length(); i++ ) {
                        dls ( i ) = new datalink_m2m;
                        dls(i)->set_connection( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), _rv(), _rvc() );
                }

                for ( int i=0; i<mpdfs.length(); i++ ) {
                        epdfs(i) = mpdfs(i)->e();
                }
        };

        double evallogcond ( const vec &val, const vec &cond ) {
                int i;
                double res = 0.0;
                for ( i = mpdfs.length() - 1;i >= 0;i-- ) {
                        /*                      if ( mpdfs(i)->_rvc().count() >0) {
                                                        mpdfs ( i )->condition ( dls ( i )->get_cond ( val,cond ) );
                                                }
                                                // add logarithms
                                                res += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );*/
                        res += mpdfs ( i )->evallogcond (
                                   dls ( i )->pushdown ( val ),
                                   dls ( i )->get_cond ( val, cond )
                               );
                }
                return res;
        }
        vec evallogcond_m(const mat &Dt, const vec &cond) {
                vec tmp(Dt.cols());
                for(int i=0;i<Dt.cols(); i++){
                        tmp(i) = evallogcond(Dt.get_col(i),cond);
                }
                return tmp;
        };
        vec evallogcond_m(const Array<vec> &Dt, const vec &cond) {
                vec tmp(Dt.length());
                for(int i=0;i<Dt.length(); i++){
                        tmp(i) = evallogcond(Dt(i),cond);
                }
                return tmp;             
        };


        //TODO smarter...
        vec samplecond ( const vec &cond ) {
            vec smp= std::numeric_limits<double>::infinity() * ones ( e()->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 = ( mpdfs.length() - 1 );i >= 0;i-- ) {
                        if ( mpdfs ( i )->dimensionc() ) {
                                mpdfs ( i )->condition ( dls ( i )->get_cond ( smp ,cond ) ); // smp is val here!!
                        }
                        smpi = epdfs ( i )->sample();
                        // copy contribution of this pdf into smp
                        dls ( i )->pushup ( smp, smpi );
                        // add ith likelihood contribution
                }
                return smp;
        }
        mat samplecond ( const vec &cond,  int N ) {
                mat Smp ( dimension(),N );
                for ( int i=0;i<N;i++ ) {Smp.set_col ( i,samplecond ( cond ) );}
                return Smp;
        }

        ~mprod() {};
        void from_setting(const Setting &set){
                Array<mpdf*> Atmp; //temporary Array
                UI::get(Atmp,set, "mpdfs", UI::compulsory);
                set_elements(Atmp,true);
        }
        
};
UIREGISTER(mprod);

class eprod: public epdf {
protected:
        Array<const epdf*> epdfs;
        Array<datalink*> dls;
public:
        eprod () : epdfs ( 0 ),dls ( 0 ) {};
        void set_parameters ( const Array<const epdf*> &epdfs0, bool named=true ) {
                epdfs=epdfs0;//.set_length ( epdfs0.length() );
                dls.set_length ( epdfs.length() );

                bool independent=true;
                if ( named ) {
                        for ( int i=0;i<epdfs.length();i++ ) {
                                independent=rv.add ( epdfs ( i )->_rv() );
                                it_assert_debug ( independent==true, "eprod:: given components are not independent." );
                        }
                        dim=rv._dsize();
                }
                else {
                        dim =0; for ( int i=0;i<epdfs.length();i++ ) {
                                dim+=epdfs ( i )->dimension();
                        }
                }
                //
                int cumdim=0;
                int dimi=0;
                int i;
                for ( i=0;i<epdfs.length();i++ ) {
                        dls ( i ) = new datalink;
                        if ( named ) {dls ( i )->set_connection ( epdfs ( i )->_rv() , rv );}
                        else {
                                dimi = epdfs ( i )->dimension();
                                dls ( i )->set_connection ( dimi, dim, linspace ( cumdim,cumdim+dimi-1 ) );
                                cumdim+=dimi;
                        }
                }
        }

        vec mean() const {
                vec tmp ( dim );
                for ( int i=0;i<epdfs.length();i++ ) {
                        vec pom = epdfs ( i )->mean();
                        dls ( i )->pushup ( tmp, pom );
                }
                return tmp;
        }
        vec variance() const {
                vec tmp ( dim ); //second moment
                for ( int i=0;i<epdfs.length();i++ ) {
                        vec pom = epdfs ( i )->mean();
                        dls ( i )->pushup ( tmp, pow ( pom,2 ) );
                }
                return tmp-pow ( mean(),2 );
        }
        vec sample() const {
                vec tmp ( dim );
                for ( int i=0;i<epdfs.length();i++ ) {
                        vec pom = epdfs ( i )->sample();
                        dls ( i )->pushup ( tmp, pom );
                }
                return tmp;
        }
        double evallog ( const vec &val ) const {
                double tmp=0;
                for ( int i=0;i<epdfs.length();i++ ) {
                        tmp+=epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );
                }
                it_assert_debug ( std::isfinite ( tmp ),"Infinite" );
                return tmp;
        }
        const epdf* operator () ( int i ) const {it_assert_debug ( i<epdfs.length(),"wrong index" );return epdfs ( i );}

        ~eprod() {for ( int i=0;i<epdfs.length();i++ ) {delete dls ( i );}}
};


class mmix : public mpdf {
protected:
        Array<mpdf*> Coms;

        shared_ptr<emix> iepdf;

public:
        mmix():iepdf(new emix()) {
            set_ep(iepdf);
        }

        void set_parameters ( const vec &w, const Array<mpdf*> &Coms ) {
                Array<epdf*> Eps ( Coms.length() );

                for ( int i = 0;i < Coms.length();i++ ) {
                        Eps(i) = Coms(i)->e();
                }

                iepdf->set_parameters(w, Eps);
        }

        void condition ( const vec &cond ) {
                for ( int i = 0;i < Coms.length();i++ ) {Coms ( i )->condition ( cond );}
        };
};

}
#endif //MX_H

---