root/bdm/estim/merger.h @ 311

/*!
  \file
  \brief Mergers for combination of pdfs
  \author Vaclav Smidl.

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

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

#ifndef MER_H
#define MER_H


#include "mixef.h"

namespace bdm
{
        using std::string;

        /*!
        @brief Function for general combination of pdfs

        Mixtures of Gaussian densities are used internally. Switching to other densities should be trivial.
        */

        class merger : public compositepdf, public epdf
        {
                protected:
                        //!Internal mixture of EF models
                        MixEF Mix;
                        //! Data link for each mpdf in mpdfs
                        Array<datalink_m2e*> dls;
                        //! Array of rvs that are not modelled by mpdfs at all (aux)
                        Array<RV> rvzs;
                        //! Data Links of rv0 mpdfs - these will be conditioned the (rv,rvc) of mpdfs
                        Array<datalink_m2e*> zdls;

                        //!Number of samples used in approximation
                        int Ns;
                        //!Number of components in a mixture
                        int Nc;
                        //!Prior on the log-normal merging model
                        double beta;
                        //! Projection to empirical density
                        eEmp eSmp;
                        //! coefficient of resampling
                        double effss_coef;

                        //! debug or not debug
                        bool DBG;
                        //! debugging file
                        it_file* dbg;
                        //! Flag if the samples are fixed or not
                        bool fix_smp;
                public:
//!Default constructor
                        merger ( const Array<mpdf*> &S ) :
                                        compositepdf ( S ), epdf ( ),
                                        Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ), dls ( n ), rvzs ( n ), zdls ( n ), eSmp()
                        {
                                RV ztmp;
                                rv = getrv ( false );
                                RV rvc; setrvc ( rv,rvc ); // Extend rv by rvc!
                                rv.add ( rvc );
                                // get dimension
                                dim = rv._dsize();

                                for ( int i=0;i<n;i++ )
                                {
                                        //Establich connection between mpdfs and merger
                                        dls ( i ) = new datalink_m2e;dls ( i )->set_connection ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), rv );
                                        // find out what is missing in each mpdf
                                        ztmp= mpdfs ( i )->_rv();
                                        ztmp.add ( mpdfs ( i )->_rvc() );
                                        rvzs ( i ) =rv.subt ( ztmp );
                                        zdls ( i ) = new datalink_m2e; zdls ( i )->set_connection ( rvzs ( i ), ztmp, rv ) ;
                                };
                                //Set Default values of parameters
                                beta=2.0;
                                Ns=100;
                                Nc=10;
                                Mix.set_method ( EM );
                                DBG = false;
                                fix_smp = false;
                        }
                        //! set debug file
                        void debug_file ( const string fname ) { if ( DBG ) delete dbg; dbg = new it_file ( fname ); if ( dbg ) DBG=true;}
//! Set internal parameters used in approximation
                        void set_parameters ( double beta0, int Ns0, int Nc0, double effss_coef0=0.5 ) {beta=beta0;
                                Ns=Ns0;
                                Nc=Nc0;
                                effss_coef=effss_coef0;
                                eSmp.set_parameters ( Ns0,false );
                        }
                        void set_grid ( Array<vec> &XYZ )
                        {
                                int dim=XYZ.length(); ivec szs ( dim );
                                for(int i=0; i<dim;i++){szs=XYZ(i).length();}
                                Ns=prod(szs);
                                eSmp.set_parameters(Ns,false);
                                Array<vec> &samples=eSmp._samples();
                                eSmp._w()=ones(Ns)/Ns;
                                                
                                //set samples
                                ivec is=zeros_i(dim);//indeces of dimensions in for cycle;
                                vec smpi(dim);
                                for(int i=0; i<Ns; i++){
                                        for(int j=0; j<dim; j++){smpi(j)=XYZ(j)(is(j)); /* jty vector*/ }
                                        samples(i)=smpi;
                                        // shift indeces
                                        for (int j=0;j<dim;j++){
                                                if (is(j)==szs(j)-1) { //j-th index is full
                                                        is(j)=0; //shift back
                                                        is(j+1)++; //increase th next dimension;
                                                        if (is(j+1)<szs(j+1)-1) break;
                                                } else {
                                                        is(j)++; break;
                                                }
                                        }
                                }
                                
                                fix_smp = true;
                        }
//!Initialize the proposal density. This function must be called before merge()!
                        void init()   ////////////// NOT FINISHED
                        {
                                Array<vec> Smps ( n );
                                //Gibbs sampling
                                for ( int i=0;i<n;i++ ) {Smps ( i ) =zeros ( 0 );}
                        }
//!Create a mixture density using known proposal
                        void merge ( const epdf* g0 );
//!Create a mixture density, make sure to call init() before the first call
                        void merge () {merge ( & ( Mix.posterior() ) );};

//! Merge log-likelihood values
                        vec lognorm_merge ( mat &lW );
//! sample from merged density
//! weight w is a
                        vec sample ( ) const { return Mix.posterior().sample();}
                        double evallog ( const vec &dt ) const
                        {
                                vec dtf=ones ( dt.length() +1 );
                                dtf.set_subvector ( 0,dt );
                                return Mix.logpred ( dtf );
                        }
                        vec mean() const
                        {
                                const Vec<double> &w = eSmp._w();
                                const Array<vec> &S = eSmp._samples();
                                vec tmp=zeros ( dim );
                                for ( int i=0; i<Ns; i++ )
                                {
                                        tmp+=w ( i ) *S ( i );
                                }
                                return tmp;
                        }
                        mat covariance() const
                        {
                                const vec &w = eSmp._w();
                                const Array<vec> &S = eSmp._samples();

                                vec mea = mean();

                                cout << sum ( w ) << "," << w*w <<endl;

                                mat Tmp=zeros ( dim, dim );
                                for ( int i=0; i<Ns; i++ )
                                {
                                        Tmp+=w ( i ) *outer_product ( S ( i ), S ( i ) );
                                }
                                return Tmp-outer_product ( mea,mea );
                        }
                        vec variance() const
                        {
                                const vec &w = eSmp._w();
                                const Array<vec> &S = eSmp._samples();

                                vec tmp=zeros ( dim );
                                for ( int i=0; i<Ns; i++ )
                                {
                                        tmp+=w ( i ) *pow ( S ( i ),2 );
                                }
                                return tmp-pow ( mean(),2 );
                        }
//! for future use
                        virtual ~merger()
                        {
                                for ( int i=0; i<n; i++ )
                                {
                                        delete dls ( i );
                                        delete zdls ( i );
                                }
                                if ( DBG ) delete dbg;
                        };

//! Access function
                        MixEF& _Mix() {return Mix;}
//! Access function
                        emix* proposal() {emix* tmp=Mix.epredictor(); tmp->set_rv(rv); return tmp;}
//! Access function
                        eEmp& _Smp() {return eSmp;}
        };

}

#endif // MER_H