Changeset 379 for bdm/estim/merger.h

Timestamp:

06/17/09 23:53:11 (15 years ago)

Author:

smidl

Message:

merger restructured

Files:

• bdm/estim/merger.h (modified) (1 diff)

bdm/estim/merger.h

@@ -19 +19 @@
 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;
-                                                }
+using std::string;
+
+//!Merging methods
+enum MERGER_METHOD {ARITHMETIC = 1, GEOMETRIC = 2, LOGNORMAL = 3};
+
+/*!
+@brief Base class for general combination of pdfs on discrete support
+
+Mixtures of Gaussian densities are used internally. Switching to other densities should be trivial.
+
+The merged pdfs are expected to be of the form:
+ \f[ f(x_i|y_i),  i=1..n \f]
+where the resulting merger is a density on \f$ \cup [x_i,y_i] \f$ .
+Note that all variables will be joined.
+
+As a result of this feature, each source must be extended to common support
+\f[ f(z_i|y_i,x_i) f(x_i|y_i) f(y_i)  i=1..n \f]
+where \f$ z_i \f$ accumulate variables that were not in the original source.
+These extensions are calculated on-the-fly.
+
+However, these operations can not be performed in general. Hence, this class merges only sources on common support, \f$ y_i={}, z_i={}, \forall i \f$.
+For merging of more general cases, use offsprings merger_mix and merger_grid.
+*/
+
+class merger_base : public compositepdf, public epdf
+{
+        protected:
+                //! Data link for each mpdf in mpdfs
+                Array<datalink_m2e*> dls;
+                //! Array of rvs that are not modelled by mpdfs at all, \f$ z_i \f$
+                Array<RV> rvzs;
+                //! Data Links for extension \f$ f(z_i|x_i,y_i) \f$
+                Array<datalink_m2e*> zdls;
+                //! number of support points
+                int Npoints;
+                //! number of sources
+                int Nsources;
+
+                //! switch of the methoh used for merging
+                MERGER_METHOD METHOD;
+                //!Prior on the log-normal merging model
+                double beta;
+
+                //! Projection to empirical density (could also be piece-wise linear)
+                eEmp eSmp;
+
+                //! debug or not debug
+                bool DBG;
+
+                //! debugging file
+                it_file* dbg_file;
+        public:
+                //! \name Constructors
+                //! @{
+
+                //!Empty constructor
+                merger_base () : compositepdf() {DBG=false;dbg_file=NULL;};
+                //!Constructor from sources
+                merger_base (const Array<mpdf*> &S) {set_sources (S);};
+                //! Function setting the main internal structures 
+                void set_sources (const Array<mpdf*> &Sources) {
+                        compositepdf::set_elements (Sources);
+                        //set sizes
+                        dls.set_size (Sources.length());
+                        rvzs.set_size (Sources.length());
+                        zdls.set_size (Sources.length());
+
+                        rv = getrv (/* checkoverlap = */ false);
+                        RV rvc; setrvc (rv, rvc);  // Extend rv by 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) ;
+                        };
+                }
+                //! set debug file
+                void set_debug_file (const string fname) { 
+                        if (DBG) delete dbg_file; 
+                        dbg_file = new it_file (fname); 
+                        if (dbg_file) DBG = true;
+                }
+                //! Set internal parameters used in approximation
+                void set_method (MERGER_METHOD MTH, double beta0=0.0) {
+                        METHOD = MTH; 
+                        beta = beta0;
+                }
+                //! Set support points from a pdf by drawing N samples
+                void set_support(const epdf &overall, int N){
+                        it_assert_debug ( rv.equal ( overall._rv() ),"Incompatible parameter overall!" );
+                        eSmp.set_statistics(&overall,N);
+                        Npoints=N;
+                }
+                
+                //! Destructor
+                virtual ~merger_base() {
+                        for (int i = 0; i < Nsources; i++) {
+                                delete dls (i);
+                                delete zdls (i);
+                        }
+                        if (DBG) delete dbg_file;
+                };
+                //!@}
+                
+                //! \name Mathematical operations
+                //!@{
+                
+                //!Merge given sources in given points
+                void merge () {
+                        if (eSmp._w().length() ==0) {it_error("Empty support points use set_support" );}
+                        //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);
                                         }
                                 }
-                                
-                                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
+
+                                vec wtmp = exp (merge_points (lW));
+                                //renormalize
+                                eSmp._w() = wtmp / sum (wtmp);
+                        } else {
+                                it_error("Sources are not compatible - use merger_mix");
+                        }
+                        ;
+                };
+
+
+                //! Merge log-likelihood values in points using method specified by parameter METHOD
+                vec merge_points (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;
-                        };
-
+                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);
+                }
+                //!@}
+
+                //! \name Access to attributes
+                //! @{
+
+                //! Access function
+                eEmp& _Smp() {return eSmp;}
+                
+                //!@}
+};
+
+class merger_mix : public merger_base
+{
+        protected:
+                //!Internal mixture of EF models
+                MixEF Mix;
+                //!Number of components in a mixture
+                int Ncoms;
+                //! coefficient of resampling
+                double effss_coef;
+
+        public:
+                //!\name Constructors
+                //!@{
+                merger_mix () {};
+                merger_mix (const Array<mpdf*> &S) {set_sources(S);};
+                //! Set sources and prepare all internal structures
+                void set_sources (const Array<mpdf*> &S) {
+                        merger_base::set_sources(S);
+                        Nsources = S.length();
+                }
+                //! Set internal parameters used in approximation
+                void set_parameters (int Ncoms0=10, double effss_coef0 = 0.5) {
+                        Ncoms = Ncoms0;
+                        effss_coef = effss_coef0;
+                }
+                //!@}
+                
+                //! \name Mathematical operations
+                //!@{
+                
+                //!Merge values using mixture approximation
+                void merge ();
+
+                //! sample from the approximating mixture
+                vec sample () const { return Mix.posterior().sample();}
+                //! loglikelihood computed on mixture models
+                double evallog (const vec &dt) const {
+                        vec dtf = ones (dt.length() + 1);
+                        dtf.set_subvector (0, dt);
+                        return Mix.logpred (dtf);
+                }
+                //!@}
+
+                //!\name Access functions
+                //!@{
 //! 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;}
-        };
+                MixEF& _Mix() {return Mix;}
+                //! Access function
+                emix* proposal() {emix* tmp = Mix.epredictor(); tmp->set_rv (rv); return tmp;}
+        //! @}
+                
+};
 
 }