Changeset 477 for library/bdm/stat/merger.h

Timestamp:

08/05/09 14:40:03 (15 years ago)

Author:

mido

Message:

panove, vite, jak jsem peclivej na upravu kodu.. snad se vam bude libit:) konfigurace je v souboru /system/astylerc

Files:

• library/bdm/stat/merger.h (modified) (2 diffs)

library/bdm/stat/merger.h

@@ -17 +17 @@
 #include "../estim/mixtures.h"
 
-namespace bdm
-{
+namespace bdm {
 using std::string;
 
@@ -43 +42 @@
 */
 
-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;
-                //! Default for METHOD
-                static const MERGER_METHOD DFLT_METHOD;
-                
-                //!Prior on the log-normal merging model
-                double beta;
-                //! default for beta
-                static const double DFLT_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, bool own=false);
-
-                //! Function setting the main internal structures
-                void set_sources (const Array<mpdf*> &Sources, bool own) {
-                        compositepdf::set_elements (Sources,own);
-                        Nsources=mpdfs.length();
-                        //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) ;
-                        };
-                }
-                //! Rectangular support  each vector of XYZ specifies (begining-end) interval for each dimension. Same number of points, \c dimsize, in each dimension. 
-                void set_support (const Array<vec> &XYZ, const int dimsize) {
-                        set_support(XYZ,dimsize*ones_i(XYZ.length()));
-                }
-                //! Rectangular support  each vector of XYZ specifies (begining-end) interval for each dimension. \c gridsize specifies number of points is each dimension.
-                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);
+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;
+        //! Default for METHOD
+        static const MERGER_METHOD DFLT_METHOD;
+
+        //!Prior on the log-normal merging model
+        double beta;
+        //! default for beta
+        static const double DFLT_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, bool own = false );
+
+        //! Function setting the main internal structures
+        void set_sources ( const Array<mpdf*> &Sources, bool own ) {
+                compositepdf::set_elements ( Sources, own );
+                Nsources = mpdfs.length();
+                //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 ) ;
+                };
+        }
+        //! Rectangular support  each vector of XYZ specifies (begining-end) interval for each dimension. Same number of points, \c dimsize, in each dimension.
+        void set_support ( const Array<vec> &XYZ, const int dimsize ) {
+                set_support ( XYZ, dimsize*ones_i ( XYZ.length() ) );
+        }
+        //! Rectangular support  each vector of XYZ specifies (begining-end) interval for each dimension. \c gridsize specifies number of points is each dimension.
+        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;
                                 }
                         }
                 }
-                //! 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=DFLT_METHOD, double beta0 = DFLT_beta) {
-                        METHOD = MTH;
-                        beta = beta0;
-                }
-                //! Set support points from a pdf by drawing N samples
-                void set_support (const epdf &overall, int N) {
-                        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
-                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");
-                        }
-                };
-
-
-                //! 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 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;}
-
-                //! load from setting
-                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);
+        }
+        //! 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 = DFLT_METHOD, double beta0 = DFLT_beta ) {
+                METHOD = MTH;
+                beta = beta0;
+        }
+        //! Set support points from a pdf by drawing N samples
+        void set_support ( const epdf &overall, int N ) {
+                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
+        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 );
                                 }
                         }
-                        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");
-                }
-                //!@}
+
+                        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" );
+                }
+        };
+
+
+        //! 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 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;
+        }
+
+        //! load from setting
+        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);
-
-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 [0,1]
-                double effss_coef;
-                //! stop after niter iterations
-                int stop_niter;
-                
-                //! default value for Ncoms
-                static const int DFLT_Ncoms;
-                //! default value for efss_coef;
-                static const double DFLT_effss_coef;
-
-        public:
-                //!\name Constructors
-                //!@{
-                merger_mix () {};
-                merger_mix (const Array<mpdf*> &S,bool own=false) {set_sources (S,own);};
-                //! Set sources and prepare all internal structures
-                void set_sources (const Array<mpdf*> &S, bool own) {
-                        merger_base::set_sources (S,own);
-                        Nsources = S.length();
-                }
-                //! Set internal parameters used in approximation
-                void set_parameters (int Ncoms0 = DFLT_Ncoms, double effss_coef0 = DFLT_effss_coef) {
-                        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
-                //!@{
+UIREGISTER ( merger_base );
+
+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 [0,1]
+        double effss_coef;
+        //! stop after niter iterations
+        int stop_niter;
+
+        //! default value for Ncoms
+        static const int DFLT_Ncoms;
+        //! default value for efss_coef;
+        static const double DFLT_effss_coef;
+
+public:
+        //!\name Constructors
+        //!@{
+        merger_mix () {};
+        merger_mix ( const Array<mpdf*> &S, bool own = false ) {
+                set_sources ( S, own );
+        };
+        //! Set sources and prepare all internal structures
+        void set_sources ( const Array<mpdf*> &S, bool own ) {
+                merger_base::set_sources ( S, own );
+                Nsources = S.length();
+        }
+        //! Set internal parameters used in approximation
+        void set_parameters ( int Ncoms0 = DFLT_Ncoms, double effss_coef0 = DFLT_effss_coef ) {
+                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;}
-                //! from_settings
-                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);
-                }
-                
-                //! @}
+        MixEF& _Mix() {
+                return Mix;
+        }
+        //! Access function
+        emix* proposal() {
+                emix* tmp = Mix.epredictor();
+                tmp->set_rv ( rv );
+                return tmp;
+        }
+        //! from_settings
+        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);
+UIREGISTER ( merger_mix );
 
 }