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

Timestamp:

06/10/10 21:54:57 (14 years ago)

Author:

smidl

Message:

Changes in merger + change in loading ARX

Files:

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

library/bdm/stat/merger.cpp

@@ -5 +5 @@
 namespace bdm {
 
-merger_base::merger_base ( const Array<shared_ptr<pdf> > &S ) :
+MergerDiscrete::MergerDiscrete ( const Array<shared_ptr<pdf> > &S ) :
     Npoints ( 0 ), DBG ( false ), dbg_file ( 0 ) {
     set_sources ( S );
@@ -11 +11 @@
 
 
-void merger_base::set_sources ( const Array<shared_ptr<pdf> > &Sources ) {
-    pdfs = Sources;
-    Nsources = pdfs.length();
+void MergerDiscrete::set_sources ( const Array<shared_ptr<pdf> > &PartSources ) {
+    part_sources = PartSources;
+    Nsources = part_sources.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 );
+    sources.set_size ( Nsources);
+        dls.set_size ( Nsources );
+        rvzs.set_size ( Nsources );
+        zdls.set_size ( Nsources );
+
+    RV rv = get_composite_rv ( part_sources, /* checkoverlap = */ false );
 
     RV rvc;
     // Extend rv by rvc!
-    for ( int i = 0; i < pdfs.length(); i++ ) {
-        RV rvx = pdfs ( i )->_rvc().subt ( rv );
+    for ( int i = 0; i < sources.length(); i++ ) {
+        RV rvx = part_sources ( i )->_rvc().subt ( rv );
         rvc.add ( rvx ); // add rv to common rvc
     }
@@ -30 +31 @@
     // 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++ ) {
+    for ( int i = 0; i < part_sources.length(); i++ ) {         
         //Establich connection between pdfs and merger
         dls ( i ) = new datalink_m2e;
-        dls ( i )->set_connection ( pdfs ( i )->_rv(), pdfs ( i )->_rvc(), rv );
+                dls ( i )->set_connection ( part_sources ( i )->_rv(), part_sources( i )->_rvc(), rv );
 
         // find out what is missing in each pdf
-        xytmp = pdfs ( i )->_rv();
-        xytmp.add ( pdfs ( i )->_rvc() );
+        xytmp = part_sources ( i )->_rv();
+                xytmp.add ( part_sources ( i )->_rvc() );
         // z_i = common_rv-xy
         rvzs ( i ) = rv.subt ( xytmp );
@@ -49 +48 @@
         zdls ( i )->set_connection ( rvzs ( i ), xytmp, rv ) ;
     };
-}
-
-void merger_base::set_support ( rectangular_support &Sup ) {
+        
+        // 
+        merger().set_rv(rv);
+}
+
+void MergerDiscrete::set_support ( rectangular_support &Sup ) {
     Npoints = Sup.points();
     eSmp.set_parameters ( Npoints, false );
@@ -61 +63 @@
         samples ( j ) = Sup.next_vec();
     }
-}
-
-void merger_base::merge () {
+    eSmp.validate();
+}
+
+void MergerDiscrete::merge () {
     validate();
 
     //check if sources overlap:
     bool OK = true;
-    for ( int i = 0; i < pdfs.length(); i++ ) {
+    for ( int i = 0; i < part_sources.length(); i++ ) {
         OK &= ( rvzs ( i )._dsize() == 0 ); // z_i is empty
-        OK &= ( pdfs ( i )->_rvc()._dsize() == 0 ); // y_i is empty
+        OK &= ( part_sources ( i )->_rvc()._dsize() == 0 ); // y_i is empty
     }
 
     if ( OK ) {
-        mat lW = zeros ( pdfs.length(), eSmp._w().length() );
+        mat lW = zeros ( part_sources.length(), eSmp._w().length() );
 
         vec emptyvec ( 0 );
-        for ( int i = 0; i < pdfs.length(); i++ ) {
+        for ( int i = 0; i < part_sources.length(); i++ ) {
             for ( int j = 0; j < eSmp._w().length(); j++ ) {
-                lW ( i, j ) = pdfs ( i )->evallogcond ( eSmp._samples() ( j ), emptyvec );
+                lW ( i, j ) = part_sources ( i )->evallogcond ( eSmp._samples() ( j ), emptyvec );
             }
         }
@@ -92 +95 @@
 }
 
-vec merger_base::merge_points ( mat &lW ) {
+vec MergerDiscrete::merge_points ( mat &lW ) {
     int nu = lW.rows();
     vec result;
@@ -144 +147 @@
 }
 
-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++ ) {
-        vec tmp=S ( i )-mea; //inefficient but numerically stable
-        Tmp += w ( i ) * outer_product (tmp , tmp );
-    }
-    return Tmp;
-}
-
-vec merger_base::variance() const {
-    return eSmp.variance();
-}
-
-
-void merger_base::from_setting ( const Setting& set ) {
+void MergerDiscrete::from_setting ( const Setting& set ) {
     // get support
     // find which method to use
-    epdf::from_setting (set);
     string meth_str;
     UI::get( meth_str, set, "method", UI::compulsory );
@@ -197 +168 @@
 }
 
-void merger_base::to_setting  (Setting  &set) const {
-    epdf::to_setting(set);
-
+void MergerDiscrete::to_setting  (Setting  &set) const {
     UI::save( METHOD, set, "method");
 
@@ -209 +178 @@
 }
 
-void merger_base::validate() {
+void MergerDiscrete::validate() {
 //        bdm_assert ( eSmp._w().length() > 0, "Empty support, use set_support()." );
 //        bdm_assert ( dim == eSmp._samples() ( 0 ).length(), "Support points and rv are not compatible!" );
-    epdf::validate();
-    bdm_assert ( isnamed(), "mergers must be named" );
-}
-
-// DEFAULTS FOR MERGER_BASE
-const MERGER_METHOD merger_base::DFLT_METHOD = LOGNORMAL;
-const double merger_base::DFLT_beta = 1.2;
+    MergerBase::validate();
+    //bdm_assert ( merger().isnamed(), "mergers must be named" );
+}
 
 void merger_mix::merge ( ) {
+        int dim = eSmp.dimension();
     if(Npoints<1) {
         set_support(enorm<fsqmat>(zeros(dim), eye(dim)), 1000);
@@ -231 +197 @@
     vec &w = eSmp._w(); //aux
 
-    mat Smp_ex = ones ( dim + 1, Npoints ); // Extended samples for the ARX model - the last row is ones
+    mat Smp_ex = ones ( dim, Npoints ); // Extended samples for the ARX model
     for ( int i = 0; i < Npoints; i++ ) {
         set_col_part ( Smp_ex, i, Smp ( i ) );
@@ -242 +208 @@
     vec lw_mix ( Npoints );        // weights of the approximating mixture
     vec lw ( Npoints );            // tmp
-    mat lW = zeros ( Nsources, Npoints ); // array of weights of all sources
+    mat lW = zeros ( Nsources, Npoints ); // array of weights of all part_sources
     vec vec0 ( 0 );
 
@@ -263 +229 @@
     char dbg_str[100];
 
-    emix* Mpred = Mix.epredictor ( );
+        emix *Mpred = Mix.epredictor ( );
     vec Mix_pdf ( Npoints );
     while ( !converged ) {
@@ -272 +238 @@
         delete Mpred;
         Mpred = Mix.epredictor ( ); // Allocation => must be deleted at the end!!
-        Mpred->set_rv ( rv ); //the predictor predicts rv of this merger
+        Mpred->set_rv ( merger()._rv() ); //the predictor predicts rv of this merger
 
         // This will be active only later in iterations!!!
@@ -299 +265 @@
             *dbg_file << Name ( dbg_str ) << Mpred->variance();
             sprintf ( dbg_str, "Mpred_cov%d", niter );
-            *dbg_file << Name ( dbg_str ) << covariance();
+            *dbg_file << Name ( dbg_str ) << Mpred->covariance();
 
 
@@ -313 +279 @@
         }
         //Importace weighting
-        for ( int i = 0; i < pdfs.length(); i++ ) {
+        for ( int i = 0; i < part_sources.length(); i++ ) {
             lw_src = 0.0;
             //======== Same RVs ===========
             //Split according to dependency in rvs
-            if ( pdfs ( i )->dimension() == dim ) {
+            if ( part_sources ( i )->dimension() == dim ) {
                 // no need for conditioning or marginalization
-                lw_src = pdfs ( i )->evallogcond_mat ( Smp , vec ( 0 ) );
+                lw_src = part_sources ( i )->evallogcond_mat ( Smp , vec ( 0 ) );
             } else {
                 // compute likelihood of marginal on the conditional variable
-                if ( pdfs ( i )->dimensionc() > 0 ) {
+                if ( part_sources ( i )->dimensionc() > 0 ) {
                     // Make marginal on rvc_i
-                    shared_ptr<epdf> tmp_marg = Mpred->marginal ( pdfs ( i )->_rvc() );
+                    shared_ptr<epdf> tmp_marg = Mpred->marginal ( part_sources ( i )->_rvc() );
                     //compute vector of lw_src
                     for ( int k = 0; k < Npoints; k++ ) {
@@ -336 +302 @@
                 }
                 // Compute likelihood of the missing variable
-                if ( dim > ( pdfs ( i )->dimension() + pdfs ( i )->dimensionc() ) ) {
+                if ( dim > ( part_sources ( i )->dimension() + part_sources ( i )->dimensionc() ) ) {
                     ///////////////
                     // There are variales unknown to pdfs(i) : rvzs
@@ -355 +321 @@
                 // Compute likelihood of the partial source
                 for ( int k = 0; k < Npoints; k++ ) {
-                    lw_src ( k ) += pdfs ( i )->evallogcond ( dls ( i )->pushdown ( Smp ( k ) ),
+                    lw_src ( k ) += part_sources ( i )->evallogcond ( dls ( i )->pushdown ( Smp ( k ) ),
                                     dls ( i )->get_cond ( Smp ( k ) ) );
                 }
@@ -363 +329 @@
             lW.set_row ( i, lw_src ); // do not divide by mix
         }
-        lw = merger_base::merge_points ( lW ); //merge
+        lw = MergerDiscrete::merge_points ( lW ); //merge
 
         //Importance weighting
@@ -391 +357 @@
     }
     delete Mpred;
-//        cout << endl;
-
 }
 
 void merger_mix::from_setting ( const Setting& set ) {
-    merger_base::from_setting ( set );
+    MergerDiscrete::from_setting ( set );
     Ncoms=DFLT_Ncoms;
     UI::get( Ncoms, set, "ncoms", UI::optional );
@@ -406 +370 @@
 
 void merger_mix::to_setting  (Setting  &set) const  {
-    merger_base::to_setting(set);
+    MergerDiscrete::to_setting(set);
     UI::save( Ncoms, set, "ncoms");
     UI::save (effss_coef , set,  "effss_coef");
@@ -413 +377 @@
 
 void merger_mix::validate() {
-    merger_base::validate();
+    MergerDiscrete::validate();
     bdm_assert(Ncoms>0,"Ncoms too small");
 }