merger.h

Go to the documentation of this file.

#ifndef MERGER_H
#define MERGER_H


#include "../estim/mixtures.h"

namespace bdm {
using std::string;

enum MERGER_METHOD {ARITHMETIC = 1, GEOMETRIC = 2, LOGNORMAL = 3};

class merger_base : public epdf {
protected:
        Array<shared_ptr<mpdf> > mpdfs;

        Array<datalink_m2e*> dls;

        Array<RV> rvzs;

        Array<datalink_m2e*> zdls;

        int Npoints;

        int Nsources;

        MERGER_METHOD METHOD;
        static const MERGER_METHOD DFLT_METHOD;

        double beta;
        static const double DFLT_beta;

        eEmp eSmp;

        bool DBG;

        it_file* dbg_file;
public:

        merger_base () : Npoints(0), Nsources(0), DBG(false), dbg_file(0) {
        }

        merger_base ( const Array<shared_ptr<mpdf> > &S );

        void set_sources ( const Array<shared_ptr<mpdf> > &Sources ) {
                mpdfs = Sources;
                Nsources = mpdfs.length();
                //set sizes
                dls.set_size ( Sources.length() );
                rvzs.set_size ( Sources.length() );
                zdls.set_size ( Sources.length() );

                rv = get_composite_rv ( mpdfs, /* checkoverlap = */ false );

                RV rvc;
                // Extend rv by rvc!
                for ( int i = 0; i < mpdfs.length(); i++ ) {
                        RV rvx = mpdfs ( i )->_rvc().subt ( rv );
                        rvc.add ( rvx ); // add rv to common rvc
                }

                // join rv and rvc - see descriprion
                rv.add ( rvc );
                // get dimension
                dim = rv._dsize();

                // create links between sources and common rv
                RV xytmp;
                for ( int i = 0; i < mpdfs.length(); 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
                        xytmp = mpdfs ( i )->_rv();
                        xytmp.add ( mpdfs ( i )->_rvc() );
                        // z_i = common_rv-xy
                        rvzs ( i ) = rv.subt ( xytmp );
                        //establish connection between extension (z_i|x,y)s and common rv
                        zdls ( i ) = new datalink_m2e;
                        zdls ( i )->set_connection ( rvzs ( i ), xytmp, rv ) ;
                };
        }
        void set_support ( const Array<vec> &XYZ, const int dimsize ) {
                set_support ( XYZ, dimsize*ones_i ( XYZ.length() ) );
        }
        void set_support ( const Array<vec> &XYZ, const ivec &gridsize ) {
                int dim = XYZ.length();  //check with internal dim!!
                Npoints = prod ( gridsize );
                eSmp.set_parameters ( Npoints, false );
                Array<vec> &samples = eSmp._samples();
                eSmp._w() = ones ( Npoints ) / Npoints; //unifrom size of bins
                //set samples
                ivec ind = zeros_i ( dim );    //indeces of dimensions in for cycle;
                vec smpi ( dim );          // ith sample
                vec steps = zeros ( dim );        // ith sample
                // first corner
                for ( int j = 0; j < dim; j++ ) {
                        smpi ( j ) = XYZ ( j ) ( 0 ); /* beginning of the interval*/
                        it_assert ( gridsize ( j ) != 0.0, "Zeros in gridsize!" );
                        steps ( j ) = ( XYZ ( j ) ( 1 ) - smpi ( j ) ) / gridsize ( j );
                }
                // fill samples
                for ( int i = 0; i < Npoints; i++ ) {
                        // copy
                        samples ( i ) = smpi;
                        // go through all dimensions
                        for ( int j = 0; j < dim; j++ ) {
                                if ( ind ( j ) == gridsize ( j ) - 1 ) { //j-th index is full
                                        ind ( j ) = 0; //shift back
                                        smpi ( j ) = XYZ ( j ) ( 0 );

                                        if ( i < Npoints - 1 ) {
                                                ind ( j + 1 ) ++; //increase the next dimension;
                                                smpi ( j + 1 ) += steps ( j + 1 );
                                                break;
                                        }

                                } else {
                                        ind ( j ) ++;
                                        smpi ( j ) += steps ( j );
                                        break;
                                }
                        }
                }
        }
        void set_debug_file ( const string fname ) {
                if ( DBG ) delete dbg_file;
                dbg_file = new it_file ( fname );
                if ( dbg_file ) DBG = true;
        }
        void set_method ( MERGER_METHOD MTH = DFLT_METHOD, double beta0 = DFLT_beta ) {
                METHOD = MTH;
                beta = beta0;
        }
        void set_support ( const epdf &overall, int N ) {
                eSmp.set_statistics ( overall, N );
                Npoints = N;
        }

        virtual ~merger_base() {
                for ( int i = 0; i < Nsources; i++ ) {
                        delete dls ( i );
                        delete zdls ( i );
                }
                if ( DBG ) delete dbg_file;
        };


        virtual void merge () {
                validate();

                //check if sources overlap:
                bool OK = true;
                for ( int i = 0; i < mpdfs.length(); i++ ) {
                        OK &= ( rvzs ( i )._dsize() == 0 ); // z_i is empty
                        OK &= ( mpdfs ( i )->_rvc()._dsize() == 0 ); // y_i is empty
                }

                if ( OK ) {
                        mat lW = zeros ( mpdfs.length(), eSmp._w().length() );

                        vec emptyvec ( 0 );
                        for ( int i = 0; i < mpdfs.length(); i++ ) {
                                for ( int j = 0; j < eSmp._w().length(); j++ ) {
                                        lW ( i, j ) = mpdfs ( i )->evallogcond ( eSmp._samples() ( j ), emptyvec );
                                }
                        }

                        vec w_nn = merge_points ( lW );
                        vec wtmp = exp ( w_nn - max ( w_nn ) );
                        //renormalize
                        eSmp._w() = wtmp / sum ( wtmp );
                } else {
                        it_error ( "Sources are not compatible - use merger_mix" );
                }
        };


        vec merge_points ( mat &lW );


        vec mean() const {
                const Vec<double> &w = eSmp._w();
                const Array<vec> &S = eSmp._samples();
                vec tmp = zeros ( dim );
                for ( int i = 0; i < Npoints; 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 < Npoints; 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 < Nsources; i++ ) {
                        tmp += w ( i ) * pow ( S ( i ), 2 );
                }
                return tmp - pow ( mean(), 2 );
        }


        eEmp& _Smp() {
                return eSmp;
        }

        void from_setting ( const Setting& set ) {
                // get support
                // find which method to use
                string meth_str;
                UI::get<string> ( meth_str, set, "method", UI::compulsory );
                if ( !strcmp ( meth_str.c_str(), "arithmetic" ) )
                        set_method ( ARITHMETIC );
                else {
                        if ( !strcmp ( meth_str.c_str(), "geometric" ) )
                                set_method ( GEOMETRIC );
                        else if ( !strcmp ( meth_str.c_str(), "lognormal" ) ) {
                                set_method ( LOGNORMAL );
                                set.lookupValue ( "beta", beta );
                        }
                }
                string dbg_file;
                if ( UI::get ( dbg_file, set, "dbg_file" ) )
                        set_debug_file ( dbg_file );
                //validate() - not used
        }

        void validate() {
                it_assert ( eSmp._w().length() > 0, "Empty support, use set_support()." );
                it_assert ( dim == eSmp._samples() ( 0 ).length(), "Support points and rv are not compatible!" );
                it_assert ( isnamed(), "mergers must be named" );
        }
};
UIREGISTER ( merger_base );
SHAREDPTR ( merger_base );

class merger_mix : public merger_base {
protected:
        MixEF Mix;
        int Ncoms;
        double effss_coef;
        int stop_niter;

        static const int DFLT_Ncoms;
        static const double DFLT_effss_coef;

public:
        merger_mix ():Ncoms(0), effss_coef(0), stop_niter(0) { }

        merger_mix ( const Array<shared_ptr<mpdf> > &S ):
                Ncoms(0), effss_coef(0), stop_niter(0) {
                set_sources ( S );
        }

        void set_sources ( const Array<shared_ptr<mpdf> > &S ) {
                merger_base::set_sources ( S );
                Nsources = S.length();
        }

        void set_parameters ( int Ncoms0 = DFLT_Ncoms, double effss_coef0 = DFLT_effss_coef ) {
                Ncoms = Ncoms0;
                effss_coef = effss_coef0;
        }


        void merge ();

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

        MixEF& _Mix() {
                return Mix;
        }
        emix* proposal() {
                emix* tmp = Mix.epredictor();
                tmp->set_rv ( rv );
                return tmp;
        }
        void from_setting ( const Setting& set ) {
                merger_base::from_setting ( set );
                set.lookupValue ( "ncoms", Ncoms );
                set.lookupValue ( "effss_coef", effss_coef );
                set.lookupValue ( "stop_niter", stop_niter );
        }


};
UIREGISTER ( merger_mix );
SHAREDPTR ( merger_mix );

}

#endif // MER_H

---