Changeset 739 for library/bdm/stat

Timestamp:

11/25/09 18:02:21 (15 years ago)

Author:

mido

Message:

the rest of h to cpp movements, with exception of from_setting and validate to avoid conflicts with Sarka

Location:

library/bdm/stat

Files:

• emix.cpp (modified) (4 diffs)
• emix.h (modified) (3 diffs)
• exp_family.cpp (modified) (2 diffs)
• exp_family.h (modified) (4 diffs)
• merger.cpp (modified) (2 diffs)
• merger.h (modified) (3 diffs)

library/bdm/stat/emix.cpp

@@ -34 +34 @@
 
         return Coms ( i )->sample();
+}
+
+vec emix::mean() const {
+        int i;
+        vec mu = zeros ( dim );
+        for ( i = 0; i < w.length(); i++ ) {
+                mu += w ( i ) * Coms ( i )->mean();
+        }
+        return mu;
+}
+
+vec emix::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 emix::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 );
+        bdm_assert_debug ( std::isfinite ( tmp ), "Infinite" );
+        return tmp;
+}
+
+vec emix::evallog_mat ( const mat &Val ) const {
+        vec x = zeros ( Val.cols() );
+        for ( int i = 0; i < w.length(); i++ ) {
+                x += w ( i ) * exp ( Coms ( i )->evallog_mat ( Val ) );
+        }
+        return log ( x );
+};
+
+mat emix::evallog_coms ( 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_mat ( Val ) ) );
+        }
+        return X;
 }
 
@@ -203 +252 @@
 };
 
+double mprod::evallogcond ( const vec &val, const vec &cond ) {
+        int i;
+        double res = 0.0;
+        for ( i = pdfs.length() - 1; i >= 0; i-- ) {
+                /*                      if ( pdfs(i)->_rvc().count() >0) {
+                                                pdfs ( i )->condition ( dls ( i )->get_cond ( val,cond ) );
+                                        }
+                                        // add logarithms
+                                        res += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );*/
+                res += pdfs ( i )->evallogcond (
+                           dls ( i )->pushdown ( val ),
+                           dls ( i )->get_cond ( val, cond )
+                       );
+        }
+        return res;
+}
+
+vec mprod::evallogcond_mat ( 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 mprod::evallogcond_mat ( 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;
+}
+
 void mprod::set_elements ( const Array<shared_ptr<pdf> > &mFacs ) {
         pdfs = mFacs;
@@ -236 +318 @@
 
         return Coms ( i )->samplecond ( cond );
+}
+
+vec eprod::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 eprod::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 eprod::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 eprod::evallog ( const vec &val ) const {
+        double tmp = 0;
+        for ( int i = 0; i < epdfs.length(); i++ ) {
+                tmp += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );
+        }
+        bdm_assert_debug ( std::isfinite ( tmp ), "Infinite" );
+        return tmp;
 }
 
@@ -267 +383 @@
 //              }
 //      };
+

library/bdm/stat/emix.h

@@ -129 +129 @@
 
         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 );
-                bdm_assert_debug ( std::isfinite ( tmp ), "Infinite" );
-                return tmp;
-        };
-
-        vec evallog_mat ( const mat &Val ) const {
-                vec x = zeros ( Val.cols() );
-                for ( int i = 0; i < w.length(); i++ ) {
-                        x += w ( i ) * exp ( Coms ( i )->evallog_mat ( Val ) );
-                }
-                return log ( x );
-        };
+
+        vec mean() const;
+
+        vec variance() const;
+
+        double evallog ( const vec &val ) const;
+
+        vec evallog_mat ( const mat &Val ) const;
 
         //! Auxiliary function that returns pdflog for each component
-        mat evallog_coms ( 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_mat ( Val ) ) );
-                }
-                return X;
-        };
+        mat evallog_coms ( const mat &Val ) const;
 
         shared_ptr<epdf> marginal ( const RV &rv ) const;
@@ -334 +298 @@
         void set_elements ( const Array<shared_ptr<pdf> > &mFacs );
 
-        double evallogcond ( const vec &val, const vec &cond ) {
-                int i;
-                double res = 0.0;
-                for ( i = pdfs.length() - 1; i >= 0; i-- ) {
-                        /*                      if ( pdfs(i)->_rvc().count() >0) {
-                                                        pdfs ( i )->condition ( dls ( i )->get_cond ( val,cond ) );
-                                                }
-                                                // add logarithms
-                                                res += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );*/
-                        res += pdfs ( i )->evallogcond (
-                                   dls ( i )->pushdown ( val ),
-                                   dls ( i )->get_cond ( val, cond )
-                               );
-                }
-                return res;
-        }
-        vec evallogcond_mat ( 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_mat ( 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;
-        };
-
+        double evallogcond ( const vec &val, const vec &cond );
+
+        vec evallogcond_mat ( const mat &Dt, const vec &cond );
+
+        vec evallogcond_mat ( const Array<vec> &Dt, const vec &cond );
 
         //TODO smarter...
@@ -441 +379 @@
         }
 
-        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 ) );
-                }
-                bdm_assert_debug ( std::isfinite ( tmp ), "Infinite" );
-                return tmp;
-        }
+        vec mean() const;
+
+        vec variance() const;
+
+        vec sample() const;
+
+        double evallog ( const vec &val ) const;
+
         //!access function
         const epdf* operator () ( int i ) const {

library/bdm/stat/exp_family.cpp

@@ -279 +279 @@
         M = iLsub * Lpsi;
         R = ldR.to_mat()  ;
+}
+
+void egiw::log_register ( bdm::logger& L, const string& prefix ) {
+        if ( log_level == 3 ) {
+                root::log_register ( L, prefix );
+                logrec->ids.set_length ( 2 );
+                int th_dim = dimension() - dimx * ( dimx + 1 ) / 2;
+                logrec->ids ( 0 ) = L.add_vector ( RV ( "", th_dim ), prefix + logrec->L.prefix_sep() + "mean" );
+                logrec->ids ( 1 ) = L.add_vector ( RV ( "", th_dim * th_dim ), prefix + logrec->L.prefix_sep() + "variance" );
+        } else {
+                epdf::log_register ( L, prefix );
+        }
+}
+
+void egiw::log_write() const {
+        if ( log_level == 3 ) {
+                mat M;
+                ldmat Lam;
+                ldmat Vz;
+                factorize ( M, Vz, Lam );
+                logrec->L.log_vector ( logrec->ids ( 0 ), est_theta() );
+                logrec->L.log_vector ( logrec->ids ( 1 ), cvectorize ( est_theta_cov().to_mat() ) );
+        } else {
+                epdf::log_write();
+        }
+
+}
+
+void multiBM::bayes ( const vec &yt, const vec &cond ) {
+        if ( frg < 1.0 ) {
+                beta *= frg;
+                last_lognc = est.lognc();
+        }
+        beta += yt;
+        if ( evalll ) {
+                ll = est.lognc() - last_lognc;
+        }
+}
+
+double multiBM::logpred ( const vec &yt ) const {
+        eDirich pred ( est );
+        vec &beta = pred._beta();
+
+        double lll;
+        if ( frg < 1.0 ) {
+                beta *= frg;
+                lll = pred.lognc();
+        } else if ( evalll ) {
+                lll = last_lognc;
+        } else {
+                lll = pred.lognc();
+        }
+
+        beta += yt;
+        return pred.lognc() - lll;
+}
+void multiBM::flatten ( const BMEF* B ) {
+        const multiBM* E = dynamic_cast<const multiBM*> ( B );
+        // sum(beta) should be equal to sum(B.beta)
+        const vec &Eb = E->beta;//const_cast<multiBM*> ( E )->_beta();
+        beta *= ( sum ( Eb ) / sum ( beta ) );
+        if ( evalll ) {
+                last_lognc = est.lognc();
+        }
 }
 
@@ -462 +526 @@
 }
 
+void mlstudent::condition ( const vec &cond ) {
+        if ( cond.length() > 0 ) {
+                iepdf._mu() = A * cond + mu_const;
+        } else {
+                iepdf._mu() =  mu_const;
+        }
+        double zeta;
+        //ugly hack!
+        if ( ( cond.length() + 1 ) == Lambda.rows() ) {
+                zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) );
+        } else {
+                zeta = Lambda.invqform ( cond );
+        }
+        _R = Re;
+        _R *= ( 1 + zeta );// / ( nu ); << nu is in Re!!!!!!
+}
+
+void eEmp::qbounds ( vec &lb, vec &ub, double perc ) const {
+        // lb in inf so than it will be pushed below;
+        lb.set_size ( dim );
+        ub.set_size ( dim );
+        lb = std::numeric_limits<double>::infinity();
+        ub = -std::numeric_limits<double>::infinity();
+        int j;
+        for ( int i = 0; i < n; i++ ) {
+                for ( j = 0; j < dim; j++ ) {
+                        if ( samples ( i ) ( j ) < lb ( j ) ) {
+                                lb ( j ) = samples ( i ) ( j );
+                        }
+                        if ( samples ( i ) ( j ) > ub ( j ) ) {
+                                ub ( j ) = samples ( i ) ( j );
+                        }
+                }
+        }
+}
+
+
 };

library/bdm/stat/exp_family.h

@@ -323 +323 @@
                 // check sizes, rvs etc.
         }
-        void log_register ( bdm::logger& L, const string& prefix ) {
-                if ( log_level == 3 ) {
-                        root::log_register ( L, prefix );
-                        logrec->ids.set_length ( 2 );
-                        int th_dim = dimension() - dimx * ( dimx + 1 ) / 2;
-                        logrec->ids ( 0 ) = L.add_vector ( RV ( "", th_dim ), prefix + logrec->L.prefix_sep() + "mean" );
-                        logrec->ids ( 1 ) = L.add_vector ( RV ( "", th_dim * th_dim ), prefix + logrec->L.prefix_sep() + "variance" );
-                } else {
-                        epdf::log_register ( L, prefix );
-                }
-        }
-        void log_write() const {
-                if ( log_level == 3 ) {
-                        mat M;
-                        ldmat Lam;
-                        ldmat Vz;
-                        factorize ( M, Vz, Lam );
-                        logrec->L.log_vector ( logrec->ids ( 0 ), est_theta() );
-                        logrec->L.log_vector ( logrec->ids ( 1 ), cvectorize ( est_theta_cov().to_mat() ) );
-                } else {
-                        epdf::log_write();
-                }
-
-        }
+        void log_register ( bdm::logger& L, const string& prefix );
+
+        void log_write() const;
         //!@}
 };
@@ -526 +505 @@
                 beta = mB->beta;
         }
-        void bayes ( const vec &yt, const vec &cond = empty_vec ) {
-                if ( frg < 1.0 ) {
-                        beta *= frg;
-                        last_lognc = est.lognc();
-                }
-                beta += yt;
-                if ( evalll ) {
-                        ll = est.lognc() - last_lognc;
-                }
-        }
-        double logpred ( const vec &yt ) const {
-                eDirich pred ( est );
-                vec &beta = pred._beta();
-
-                double lll;
-                if ( frg < 1.0 ) {
-                        beta *= frg;
-                        lll = pred.lognc();
-                } else if ( evalll ) {
-                        lll = last_lognc;
-                } else {
-                        lll = pred.lognc();
-                }
-
-                beta += yt;
-                return pred.lognc() - lll;
-        }
-        void flatten ( const BMEF* B ) {
-                const multiBM* E = dynamic_cast<const multiBM*> ( B );
-                // sum(beta) should be equal to sum(B.beta)
-                const vec &Eb = E->beta;//const_cast<multiBM*> ( E )->_beta();
-                beta *= ( sum ( Eb ) / sum ( beta ) );
-                if ( evalll ) {
-                        last_lognc = est.lognc();
-                }
-        }
+        void bayes ( const vec &yt, const vec &cond = empty_vec );
+
+        double logpred ( const vec &yt ) const;
+
+        void flatten ( const BMEF* B );
+
         //! return correctly typed posterior (covariant return)
         const eDirich& posterior() const {
@@ -971 +920 @@
                 Lambda = Lambda0;
         }
-        void condition ( const vec &cond ) {
-                if ( cond.length() > 0 ) {
-                        iepdf._mu() = A * cond + mu_const;
-                } else {
-                        iepdf._mu() =  mu_const;
-                }
-                double zeta;
-                //ugly hack!
-                if ( ( cond.length() + 1 ) == Lambda.rows() ) {
-                        zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) );
-                } else {
-                        zeta = Lambda.invqform ( cond );
-                }
-                _R = Re;
-                _R *= ( 1 + zeta );// / ( nu ); << nu is in Re!!!!!!
-        };
+
+        void condition ( const vec &cond ); 
 
         void validate() {
@@ -1526 +1461 @@
         }
         //! For this class, qbounds are minimum and maximum value of the population!
-        void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const {
-                // lb in inf so than it will be pushed below;
-                lb.set_size ( dim );
-                ub.set_size ( dim );
-                lb = std::numeric_limits<double>::infinity();
-                ub = -std::numeric_limits<double>::infinity();
-                int j;
-                for ( int i = 0; i < n; i++ ) {
-                        for ( j = 0; j < dim; j++ ) {
-                                if ( samples ( i ) ( j ) < lb ( j ) ) {
-                                        lb ( j ) = samples ( i ) ( j );
-                                }
-                                if ( samples ( i ) ( j ) > ub ( j ) ) {
-                                        ub ( j ) = samples ( i ) ( j );
-                                }
-                        }
-                }
-        }
+        void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const;
 };
 

library/bdm/stat/merger.cpp

@@ -8 +8 @@
                 Npoints ( 0 ), DBG ( false ), dbg_file ( 0 ) {
         set_sources ( S );
+}
+
+void merger_base::set_sources ( const Array<shared_ptr<pdf> > &Sources ) {
+        pdfs = Sources;
+        Nsources = pdfs.length();
+        //set sizes
+        dls.set_size ( Sources.length() );
+        rvzs.set_size ( Sources.length() );
+        zdls.set_size ( Sources.length() );
+
+        rv = get_composite_rv ( pdfs, /* checkoverlap = */ false );
+
+        RV rvc;
+        // Extend rv by rvc!
+        for ( int i = 0; i < pdfs.length(); i++ ) {
+                RV rvx = pdfs ( 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 < pdfs.length(); i++ ) {
+                //Establich connection between pdfs and merger
+                dls ( i ) = new datalink_m2e;
+                dls ( i )->set_connection ( pdfs ( i )->_rv(), pdfs ( i )->_rvc(), rv );
+
+                // find out what is missing in each pdf
+                xytmp = pdfs ( i )->_rv();
+                xytmp.add ( pdfs ( 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 merger_base::set_support ( rectangular_support &Sup ) {
+        Npoints = Sup.points();
+        eSmp.set_parameters ( Npoints, false );
+        Array<vec> &samples = eSmp._samples();
+        eSmp._w() = ones ( Npoints ) / Npoints; //unifrom size of bins
+        //set samples
+        samples ( 0 ) = Sup.first_vec();
+        for ( int j = 1; j < Npoints; j++ ) {
+                samples ( j ) = Sup.next_vec();
+        }
+}
+
+void merger_base::merge () {
+        validate();
+
+        //check if sources overlap:
+        bool OK = true;
+        for ( int i = 0; i < pdfs.length(); i++ ) {
+                OK &= ( rvzs ( i )._dsize() == 0 ); // z_i is empty
+                OK &= ( pdfs ( i )->_rvc()._dsize() == 0 ); // y_i is empty
+        }
+
+        if ( OK ) {
+                mat lW = zeros ( pdfs.length(), eSmp._w().length() );
+
+                vec emptyvec ( 0 );
+                for ( int i = 0; i < pdfs.length(); i++ ) {
+                        for ( int j = 0; j < eSmp._w().length(); j++ ) {
+                                lW ( i, j ) = pdfs ( 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 {
+                bdm_error ( "Sources are not compatible - use merger_mix" );
+        }
 }
 
@@ -62 +143 @@
 }
 
+vec merger_base::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 merger_base::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 merger_base::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 );
+}
+
 void merger_mix::merge ( ) {
         Array<vec> &Smp = eSmp._samples(); //aux

library/bdm/stat/merger.h

@@ -93 +93 @@
 
         //! Function setting the main internal structures
-        void set_sources ( const Array<shared_ptr<pdf> > &Sources ) {
-                pdfs = Sources;
-                Nsources = pdfs.length();
-                //set sizes
-                dls.set_size ( Sources.length() );
-                rvzs.set_size ( Sources.length() );
-                zdls.set_size ( Sources.length() );
-
-                rv = get_composite_rv ( pdfs, /* checkoverlap = */ false );
-
-                RV rvc;
-                // Extend rv by rvc!
-                for ( int i = 0; i < pdfs.length(); i++ ) {
-                        RV rvx = pdfs ( 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 < pdfs.length(); i++ ) {
-                        //Establich connection between pdfs and merger
-                        dls ( i ) = new datalink_m2e;
-                        dls ( i )->set_connection ( pdfs ( i )->_rv(), pdfs ( i )->_rvc(), rv );
-
-                        // find out what is missing in each pdf
-                        xytmp = pdfs ( i )->_rv();
-                        xytmp.add ( pdfs ( 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_sources ( const Array<shared_ptr<pdf> > &Sources );
+
         //! Set support points from rectangular grid
-        void set_support ( rectangular_support &Sup ) {
-                Npoints = Sup.points();
-                eSmp.set_parameters ( Npoints, false );
-                Array<vec> &samples = eSmp._samples();
-                eSmp._w() = ones ( Npoints ) / Npoints; //unifrom size of bins
-                //set samples
-                samples ( 0 ) = Sup.first_vec();
-                for ( int j = 1; j < Npoints; j++ ) {
-                        samples ( j ) = Sup.next_vec();
-                }
-        }
+        void set_support ( rectangular_support &Sup );
+
         //! Set support points from dicrete grid
         void set_support ( discrete_support &Sup ) {
@@ -180 +134 @@
 
         //!Merge given sources in given points
-        virtual void merge () {
-                validate();
-
-                //check if sources overlap:
-                bool OK = true;
-                for ( int i = 0; i < pdfs.length(); i++ ) {
-                        OK &= ( rvzs ( i )._dsize() == 0 ); // z_i is empty
-                        OK &= ( pdfs ( i )->_rvc()._dsize() == 0 ); // y_i is empty
-                }
-
-                if ( OK ) {
-                        mat lW = zeros ( pdfs.length(), eSmp._w().length() );
-
-                        vec emptyvec ( 0 );
-                        for ( int i = 0; i < pdfs.length(); i++ ) {
-                                for ( int j = 0; j < eSmp._w().length(); j++ ) {
-                                        lW ( i, j ) = pdfs ( 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 {
-                        bdm_error ( "Sources are not compatible - use merger_mix" );
-                }
-        };
-
+        virtual void merge ();
 
         //! Merge log-likelihood values in points using method specified by parameter METHOD
@@ -216 +142 @@
         //! 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 );
-        }
+        vec mean() const;
+
+        mat covariance() const;
+
+        vec variance() const;
+
         //!@}