Changeset 739 for library/bdm/stat/merger.cpp

Timestamp:

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

Author:

mido

Message:

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

Files:

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

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