- Main Page
- Related Pages
- Modules
- Classes
- Files

- File List
- File Members

# merger.h

Go to the documentation of this file.

    00001
    00013 #ifndef MER_H
    00014 #define MER_H
    00015
    00016
    00017 #include "mixef.h"
    00018
    00019 namespace bdm
    00020 {
    00021         using std::string;
    00022
    00029         class merger : public compositepdf, public epdf
    00030         {
    00031                 protected:
    00033                         MixEF Mix;
    00035                         Array<datalink_m2e*> dls;
    00037                         Array<RV> rvzs;
    00039                         Array<datalink_m2e*> zdls;
    00040
    00042                         int Ns;
    00044                         int Nc;
    00046                         double beta;
    00048                         eEmp eSmp;
    00049
    00051                         bool DBG;
    00053                         it_file* dbg;
    00055                         bool fix_smp;
    00056                 public:
    00058                         merger ( const Array<mpdf*> &S ) :
    00059                                         compositepdf ( S ), epdf ( ),
    00060                                         Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ), dls ( n ), rvzs ( n ), zdls ( n ), eSmp()
    00061                         {
    00062                                 RV ztmp;
    00063                                 rv = getrv ( false );
    00064                                 RV rvc; setrvc ( rv,rvc ); // Extend rv by rvc!
    00065                                 rv.add ( rvc );
    00066                                 // get dimension
    00067                                 dim = rv._dsize();
    00068
    00069                                 for ( int i=0;i<n;i++ )
    00070                                 {
    00071                                         //Establich connection between mpdfs and merger
    00072                                         dls ( i ) = new datalink_m2e;dls ( i )->set_connection ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), rv );
    00073                                         // find out what is missing in each mpdf
    00074                                         ztmp= mpdfs ( i )->_rv();
    00075                                         ztmp.add ( mpdfs ( i )->_rvc() );
    00076                                         rvzs ( i ) =rv.subt ( ztmp );
    00077                                         zdls ( i ) = new datalink_m2e; zdls ( i )->set_connection ( rvzs ( i ), ztmp, rv ) ;
    00078                                 };
    00079                                 //Set Default values of parameters
    00080                                 beta=2.0;
    00081                                 Ns=100;
    00082                                 Nc=10;
    00083                                 Mix.set_method ( EM );
    00084                                 DBG = false;
    00085                                 fix_smp = false;
    00086                         }
    00088                         void debug_file ( const string fname ) { if ( DBG ) delete dbg; dbg = new it_file ( fname ); if ( dbg ) DBG=true;}
    00090                         void set_parameters ( double beta0, int Ns0, int Nc0 ) {beta=beta0;Ns=Ns0;Nc=Nc0;eSmp.set_parameters ( Ns0,false );}
    00091                         void set_grid ( Array<vec> &XYZ )
    00092                         {
    00093                                 int dim=XYZ.length(); ivec szs ( dim );
    00094                                 for(int i=0; i<dim;i++){szs=XYZ(i).length();}
    00095                                 Ns=prod(szs);
    00096                                 eSmp.set_parameters(Ns,false);
    00097                                 Array<vec> &samples=eSmp._samples();
    00098                                 eSmp._w()=ones(Ns)/Ns;
    00099
    00100                                 //set samples
    00101                                 ivec is=zeros_i(dim);//indeces of dimensions in for cycle;
    00102                                 vec smpi(dim);
    00103                                 for(int i=0; i<Ns; i++){
    00104                                         for(int j=0; j<dim; j++){smpi(j)=XYZ(j)(is(j)); /* jty vector*/ }
    00105                                         samples(i)=smpi;
    00106                                         // shift indeces
    00107                                         for (int j=0;j<dim;j++){
    00108                                                 if (is(j)==szs(j)-1) { //j-th index is full
    00109                                                         is(j)=0; //shift back
    00110                                                         is(j+1)++; //increase th next dimension;
    00111                                                         if (is(j+1)<szs(j+1)-1) break;
    00112                                                 } else {
    00113                                                         is(j)++; break;
    00114                                                 }
    00115                                         }
    00116                                 }
    00117
    00118                                 fix_smp = true;
    00119                         }
    00121                         void init()
    00122                         {
    00123                                 Array<vec> Smps ( n );
    00124                                 //Gibbs sampling
    00125                                 for ( int i=0;i<n;i++ ) {Smps ( i ) =zeros ( 0 );}
    00126                         }
    00128                         void merge ( const epdf* g0 );
    00130                         void merge () {merge ( & ( Mix.posterior() ) );};
    00131
    00133                         vec lognorm_merge ( mat &lW );
    00136                         vec sample ( ) const { return Mix.posterior().sample();}
    00137                         double evallog ( const vec &dt ) const
    00138                         {
    00139                                 vec dtf=ones ( dt.length() +1 );
    00140                                 dtf.set_subvector ( 0,dt );
    00141                                 return Mix.logpred ( dtf );
    00142                         }
    00143                         vec mean() const
    00144                         {
    00145                                 const Vec<double> &w = eSmp._w();
    00146                                 const Array<vec> &S = eSmp._samples();
    00147                                 vec tmp=zeros ( dim );
    00148                                 for ( int i=0; i<Ns; i++ )
    00149                                 {
    00150                                         tmp+=w ( i ) *S ( i );
    00151                                 }
    00152                                 return tmp;
    00153                         }
    00154                         mat covariance() const
    00155                         {
    00156                                 const vec &w = eSmp._w();
    00157                                 const Array<vec> &S = eSmp._samples();
    00158
    00159                                 vec mea = mean();
    00160
    00161                                 cout << sum ( w ) << "," << w*w <<endl;
    00162
    00163                                 mat Tmp=zeros ( dim, dim );
    00164                                 for ( int i=0; i<Ns; i++ )
    00165                                 {
    00166                                         Tmp+=w ( i ) *outer_product ( S ( i ), S ( i ) );
    00167                                 }
    00168                                 return Tmp-outer_product ( mea,mea );
    00169                         }
    00170                         vec variance() const
    00171                         {
    00172                                 const vec &w = eSmp._w();
    00173                                 const Array<vec> &S = eSmp._samples();
    00174
    00175                                 vec tmp=zeros ( dim );
    00176                                 for ( int i=0; i<Ns; i++ )
    00177                                 {
    00178                                         tmp+=w ( i ) *pow ( S ( i ),2 );
    00179                                 }
    00180                                 return tmp-pow ( mean(),2 );
    00181                         }
    00183                         virtual ~merger()
    00184                         {
    00185                                 for ( int i=0; i<n; i++ )
    00186                                 {
    00187                                         delete dls ( i );
    00188                                         delete zdls ( i );
    00189                                 }
    00190                                 if ( DBG ) delete dbg;
    00191                         };
    00192
    00194                         MixEF& _Mix() {return Mix;}
    00196                         eEmp& _Smp() {return eSmp;}
    00197         };
    00198
    00199 }
    00200
    00201 #endif // MER_H