Changeset 1079 for library

Timestamp:

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

Author:

smidl

Message:

Changes in merger + change in loading ARX

Location:

library

Files:

• bdm/base/user_info.cpp (modified) (1 diff)
• bdm/math/square_mat.h (modified) (1 diff)
• bdm/mpdm/arx_agent.h (modified) (1 diff)
• bdm/stat/exp_family.cpp (modified) (2 diffs)
• bdm/stat/exp_family.h (modified) (1 diff)
• bdm/stat/merger.cpp (modified) (21 diffs)
• bdm/stat/merger.h (modified) (12 diffs)
• tests/testsuite/merger.cfg (modified) (2 diffs)
• tests/testsuite/merger_test.cpp (modified) (3 diffs)

library/bdm/base/user_info.cpp

@@ -382 +382 @@
     if(element.exists("L")) {
         UI::from_setting(matrix.__L(), element["L"]);
-    }
+        } else { 
+                throw  UISettingException ( "UIException: ldmat must contain matrix L.", element);
+        }
     if(element.exists("D")) {
         UI::from_setting(matrix.__D(), element["D"]);
-    }
+    } else {
+                throw  UISettingException ( "UIException: ldmat must contain vector D.", element);
+        }
     matrix.validate();
 }

library/bdm/math/square_mat.h

@@ -324 +324 @@
     }
     void validate() {
+                bdm_assert(L.rows()==D.length(),"Incompatible L and D in ldmat");
         dim= L.rows();
     }

library/bdm/mpdm/arx_agent.h

@@ -86 +86 @@
             merger->merge();
             enorm<chmat> joint_pred;
-            mat Cov=merger->covariance();
+            mat Cov=merger->merger().covariance();
             if (sumsum(Cov)==0.0) {
                 bdm_error("merging failed");
             }
-            joint_pred.set_parameters(merger->mean(), Cov);
-            joint_pred.set_rv(merger->_rv());
+            joint_pred.set_parameters(merger->merger().mean(), Cov);
+            joint_pred.set_rv(merger->merger()._rv());
 
             enorm<chmat> marg_pred; // remove rvs taht should not be there

library/bdm/stat/exp_family.cpp

@@ -323 +323 @@
     }
     mat Vful;
-    if (!UI::get(V, set, "V", UI::optional)) {
-        if ( !UI::get ( Vful, set, "V", UI::optional ) ) {
+        if (set.exists("V.L")){
+                UI::get(V, set, "V", UI::optional);
+        } else {
+        if ( !UI::get ( Vful, set, "fV", UI::optional ) ) {
             vec dV;
             UI::get ( dV, set, "dV", UI::compulsory );
@@ -346 +348 @@
     nPsi = V.rows() - dimx;
     dim = dimx * ( dimx + nPsi );
-
+        if (dim<dimx) {bdm_error("Check if matrix V has correct dimension");}
+        
     if ( nu < 0 ) {
         nu = 0.1 + nPsi + 2 * dimx + 2; // +2 assures finite expected value of R

library/bdm/stat/exp_family.h

@@ -484 +484 @@
     V.L     = [...];       % L part of matrix V
     V.D     = [...];       % D part of matrix V
-    -or- V  = [...];       % full matrix V
+    -or- fV = [...];       % full matrix V
     -or- dV = [...];       % vector of diagonal of V (when V not given)
 

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");
 }

library/bdm/stat/merger.h

@@ -32 +32 @@
                 //! weights of the sources -- no restrictions
                 vec w;
-                //! get ith source
+                //! source pdfs -- conditioning allowed
                 Array<shared_ptr<epdf> > sources;
                 //! merge 
-                void merge() NOT_IMPLEMENTED_VOID;
+                virtual void merge() NOT_IMPLEMENTED_VOID;
                 //! 
                 //! check if all epdfs are on the same support
-                void validate() {
+                void validate_sources() {
                         int n=sources.length();
                         
@@ -75 +75 @@
                         vec mea=zeros(dim);
                         // go through all sources
+                        epdf * si;
                         for (int i=0; i<n; i++){
-                                sq_T Ci = sources(i)->covariance();
+                                si = dynamic_cast<epdf*>(sources(i).get());  // transform pdf into epdf
+
+                                if (!si) {bdm_error("Can't merger this type of pdf: " + sources(i)->to_string());}
+                                
+                                sq_T Ci = si->covariance();
                                 sq_T wCi = Ci; wCi*=w(i);
-                                vec mi = sources(i)->mean();
+                                vec mi = si->mean();
                                 switch (method) {
                                         case ARITHMETIC:{
@@ -133 +138 @@
 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.
+
+\TODO use sources for extensions
 */
-class merger_base : public epdf {
+class MergerDiscrete : public MergerBase {
 protected:
-    //! Elements of composition
-    Array<shared_ptr<pdf> > pdfs;
-
+        //! partial sources
+        Array<shared_ptr< pdf > > part_sources;
+        
     //! Data link for each pdf in pdfs
     Array<datalink_m2e*> dls;
@@ -177 +184 @@
 
     //! Default constructor
-    merger_base () : Npoints ( 0 ), Nsources ( 0 ), DBG ( false ), dbg_file ( 0 ) {
+    MergerDiscrete () : Npoints ( 0 ), Nsources ( 0 ), DBG ( false ), dbg_file ( 0 ) {
     }
 
     //!Constructor from sources
-    merger_base ( const Array<shared_ptr<pdf> > &S );
+    MergerDiscrete ( const Array<shared_ptr<pdf> > &S );
 
     //! Function setting the main internal structures
-    void set_sources ( const Array<shared_ptr<pdf> > &Sources );
+    void set_sources ( const Array<shared_ptr<pdf> > &PartSources );
 
     //! Set support points from rectangular grid
@@ -205 +212 @@
 
     //! Set internal parameters used in approximation
-    void set_method ( MERGER_METHOD MTH = DFLT_METHOD, double beta0 = DFLT_beta ) {
+    void set_method ( MERGER_METHOD MTH = GEOMETRIC, double beta0 = 1.2 ) {
         METHOD = MTH;
         beta = beta0;
@@ -216 +223 @@
 
     //! Destructor
-    virtual ~merger_base() {
+    virtual ~MergerDiscrete() {
         for ( int i = 0; i < Nsources; i++ ) {
             delete dls ( i );
@@ -229 +236 @@
 
     //!Merge given sources in given points
-    virtual void merge ();
+    void merge ();
 
     //! 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;
-
-    mat covariance() const;
-
-    vec variance() const;
-
-    //! Compute log-probability of argument \c val
-    virtual double evallog ( const vec &val ) const NOT_IMPLEMENTED(0);
-
-    //! Returns a sample, \f$ x \f$ from density \f$ f_x()\f$
-    virtual vec sample() const NOT_IMPLEMENTED(0);
-
-    //!@}
+  //!@}
 
     //! \name Access to attributes
@@ -258 +250 @@
         return eSmp;
     }
+
+        eEmp & merger() {return eSmp;}
 
     /*! Create object from the following structure
@@ -279 +273 @@
     //!@}
 };
-UIREGISTER ( merger_base );
-SHAREDPTR ( merger_base );
+UIREGISTER ( MergerDiscrete );
+SHAREDPTR ( MergerDiscrete );
 
 //! \brief Merger using importance sampling with mixture proposal density
-class merger_mix : public merger_base {
+class merger_mix : public MergerDiscrete {
 protected:
     //!Internal mixture of EF models
@@ -311 +305 @@
     //! Set sources and prepare all internal structures
     void set_sources ( const Array<shared_ptr<pdf> > &S ) {
-        merger_base::set_sources ( S );
+        MergerDiscrete::set_sources ( S );
         //Nsources = S.length();
     }
@@ -349 +343 @@
     emix* proposal() {
         emix* tmp = Mix.epredictor();
-        tmp->set_rv ( rv );
+        tmp->set_rv ( merger()._rv() );
         return tmp;
     }
@@ -376 +370 @@
     void validate();
 
+//      const emix &merger(){return *Mpred;}
     //! @}
 };

library/tests/testsuite/merger.cfg

@@ -12 +12 @@
   {
       class = "egiw";
-      V = ( "matrix", 2, 2, [ 20.0, 8.0, 8.0, 4.0 ] );
+      fV = ( "matrix", 2, 2, [ 20.0, 8.0, 8.0, 4.0 ] );
       nu = 4.0;
       dimx = 1.0;
@@ -29 +29 @@
 Merger : 
 {
-  class = "merger_base";
+  class = "MergerDiscrete";
   beta = 1.0;
   method = "geometric";

library/tests/testsuite/merger_test.cpp

@@ -47 +47 @@
     lW.set_row ( 1, l_f2 );
 
-    merger_base M ( A );
+    MergerDiscrete M ( A );
     enorm<fsqmat> g0;
     g0.set_rv ( x );
@@ -153 +153 @@
     UIFile in ( "merger.cfg" );
 
-    shared_ptr<merger_base> mb =
-        UI::build<merger_base> ( in, "Merger", UI::compulsory );
+    shared_ptr<MergerDiscrete> mb =
+        UI::build<MergerDiscrete> ( in, "Merger", UI::compulsory );
 
     pdf_array sources;
@@ -165 +165 @@
 
     mb->merge();
-    vec m = mb->mean();
+    vec m = mb->merger().mean();
     CHECK_EQUAL ( 2, m.size() );
 }