Changeset 190

Timestamp:

10/22/08 10:46:38 (16 years ago)

Author:

smidl

Message:

adaptation of merger for changes and creation of datalink class

Files:

• bdm/estim/merger.cpp (modified) (4 diffs)
• bdm/estim/merger.h (modified) (1 diff)
• bdm/stat/libBM.h (modified) (17 diffs)
• matlab/merger_iter_debug.m (modified) (2 diffs)
• tests/CMakeLists.txt (modified) (1 diff)
• tests/datalink_test.cpp (added)
• tests/merger_iter_test.cpp (modified) (4 diffs)

bdm/estim/merger.cpp

@@ -33 +33 @@
         vec &w = eSmp._w(); //aux
 
-        mat Smp_ex =ones ( rv.count() +1,Ns ); // Extended samples The last row is ones for the ARX model
+        mat Smp_ex =ones ( rv.count() +1,Ns ); // Extended samples for the ARX model - the last row is ones 
         for ( int i=0;i<Ns;i++ ) {      set_col_part ( Smp_ex,i,Smp ( i ) );}
 
@@ -49 +49 @@
                  V0, rv.count() *rv.count() +3.0 ); //initial guess of Mix: zero mean, large variance
 
-
-
         // ============= MAIN LOOP ==================
         bool converged=false;
@@ -62 +60 @@
                 //Re-Initialize Mixture model
                 Mix.init ( &A0, Smp_ex, Nc );
-                Mix.bayesB ( Smp_ex );
+                Mix.bayesB ( Smp_ex, ones(Ns));//w*Ns );
                 Mpred = Mix.predictor(rv); // Allocation => must be deleted at the end!!
         
@@ -152 +150 @@
                 // ==== stopping rule ===
                 niter++;
-                converged = ( niter>6);
+                converged = ( niter>20);
         }
 

bdm/estim/merger.h

@@ -41 +41 @@
                         compositepdf ( S ), epdf ( getrv ( false ) ),
                         Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) )
-        {setindices(rv); beta=2.0; Ns=100; Nc=10;}
+        {setindices(rv); beta=2.0; Ns=100; Nc=10; Mix.set_method(EM);}
         //! Set internal parameters used in approximation
         void set_parameters ( double beta0, int Ns0, int Nc0 ) {        beta=beta0;Ns=Ns0;Nc=Nc0;}

bdm/stat/libBM.h

@@ -15 +15 @@
 
 #include <itpp/itbase.h>
+#include "../itpp_ext.h"
 //#include <std>
 
@@ -20 +21 @@
 
 //! Structure of RV (used internally), i.e. expanded RVs
-class str{
+class str {
 public:
         //! vector id ids (non-unique!)
@@ -27 +28 @@
         ivec times;
         //!Default constructor
-        str(ivec ids0, ivec times0):ids(ids0),times(times0){
-                it_assert_debug(times0.length()==ids0.length(),"Incompatible input");
+        str ( ivec ids0, ivec times0 ) :ids ( ids0 ),times ( times0 ) {
+                it_assert_debug ( times0.length() ==ids0.length(),"Incompatible input" );
         };
 };
@@ -55 +56 @@
 private:
         //! auxiliary function used in constructor
-        void init (ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times );
-public:
-        //! Full constructor 
+        void init ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times );
+public:
+        //! Full constructor
         RV ( Array<std::string> in_names, ivec in_sizes, ivec in_times );
         //! Constructor with times=0
@@ -77 +78 @@
 
         //! Find indexes of self in another rv, \return ivec of the same size as self.
-        ivec findself (const RV &rv2 ) const;
+        ivec findself ( const RV &rv2 ) const;
         //! Compare if \c rv2 is identical to this \c RV
-        bool equal (const RV &rv2 ) const;
+        bool equal ( const RV &rv2 ) const;
         //! Add (concat) another variable to the current one, \return true if all rv2 were added, false if rv2 is in conflict
         bool add ( const RV &rv2 );
@@ -92 +93 @@
         //! generate \c str from rv, by expanding sizes
         str tostr() const;
-        //! when this rv is a part of bigger rv, this function returns indeces of self in the data vector of the bigger crv.  
+        //! when this rv is a part of bigger rv, this function returns indeces of self in the data vector of the bigger crv.
         //! Then, data can be copied via: data_of_this = cdata(ind);
-        ivec dataind(const RV &crv) const;
-        //! generate mutual indeces when copying data betwenn self and crv. 
+        ivec dataind ( const RV &crv ) const;
+        //! generate mutual indeces when copying data betwenn self and crv.
         //! Data are copied via: data_of_this(selfi) = data_of_rv2(rv2i)
-        void dataind(const RV &rv2, ivec &selfi, ivec &rv2i) const;
+        void dataind ( const RV &rv2, ivec &selfi, ivec &rv2i ) const;
 
         //!access function
@@ -110 +111 @@
         //!access function
         std::string name ( int at ) {return names ( at );};
-        
-        //!access function
-        void set_id ( int at, int id0 ) {ids ( at )=id0;};
-        //!access function
-        void set_size ( int at, int size0 ) {sizes ( at )=size0; tsize=sum(sizes);};
-        //!access function
-        void set_time ( int at, int time0 ) {times ( at )=time0;};
-
-        //!Assign unused ids to this rv 
+
+        //!access function
+        void set_id ( int at, int id0 ) {ids ( at ) =id0;};
+        //!access function
+        void set_size ( int at, int size0 ) {sizes ( at ) =size0; tsize=sum ( sizes );};
+        //!access function
+        void set_time ( int at, int time0 ) {times ( at ) =time0;};
+
+        //!Assign unused ids to this rv
         void newids();
 };
@@ -164 +165 @@
 //      //! Returns the required moment of the epdf
 //      virtual vec moment ( const int order = 1 );
-        
+
         //! Returns a sample, \f$x\f$ from density \f$epdf(rv)\f$
         virtual vec sample () const =0;
@@ -178 +179 @@
         virtual vec evalpdflog_m ( const mat &Val ) const {
                 vec x ( Val.cols() );
-                for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evalpdflog( Val.get_col(i) ) ;}
+                for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evalpdflog ( Val.get_col ( i ) ) ;}
                 return x;
         }
         //! Return conditional density on the given RV, the remaining rvs will be in conditioning
-        virtual mpdf* condition(const RV &rv) const  {it_warning("Not implemented"); return NULL;}
+        virtual mpdf* condition ( const RV &rv ) const  {it_warning ( "Not implemented" ); return NULL;}
         //! Return marginal density on the given RV, the remainig rvs are intergrated out
-        virtual epdf* marginal(const RV &rv) const {it_warning("Not implemented"); return NULL;}
+        virtual epdf* marginal ( const RV &rv ) const {it_warning ( "Not implemented" ); return NULL;}
 
         //! return expected value
@@ -194 +195 @@
         const RV& _rv() const {return rv;}
         //! modifier function - useful when copying epdfs
-        void _renewrv(const RV &in_rv){rv=in_rv;}
+        void _renewrv ( const RV &in_rv ) {rv=in_rv;}
         //!
 };
@@ -215 +216 @@
 //      virtual fnc moment ( const int order = 1 );
         //! Returns a sample from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv \param ll is a return value of log-likelihood of the sample.
-        virtual vec samplecond (const vec &cond, double &ll ) {this->condition ( cond );
-        vec temp= ep->sample();
-        ll=ep->evalpdflog ( temp );return temp;};
+        virtual vec samplecond ( const vec &cond, double &ll ) {
+                this->condition ( cond );
+                vec temp= ep->sample();
+                ll=ep->evalpdflog ( temp );return temp;
+        };
         //! Returns \param N samples from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv \param ll is a return value of log-likelihood of the sample.
-        virtual mat samplecond (const vec &cond, vec &ll, int N ) {
+        virtual mat samplecond ( const vec &cond, vec &ll, int N ) {
                 this->condition ( cond );
-                mat temp ( rv.count(),N ); vec smp ( rv.count() ); 
+                mat temp ( rv.count(),N ); vec smp ( rv.count() );
                 for ( int i=0;i<N;i++ ) {smp=ep->sample() ;temp.set_col ( i, smp );ll ( i ) =ep->evalpdflog ( smp );}
                 return temp;
         };
         //! Update \c ep so that it represents this mpdf conditioned on \c rvc = cond
-        virtual void condition ( const vec &cond ) {it_error("Not implemented");};
+        virtual void condition ( const vec &cond ) {it_error ( "Not implemented" );};
 
         //! Shortcut for conditioning and evaluation of the internal epdf. In some cases,  this operation can be implemented efficiently.
@@ -237 +240 @@
         mpdf ( const RV &rv0, const RV &rvc0 ) :rv ( rv0 ),rvc ( rvc0 ) {};
         //! access function
-        RV _rvc() {return rvc;}
+        RV _rvc() const {return rvc;}
         //! access function
-        RV _rv() {return rv;}
+        RV _rv() const {return rv;}
         //!access function
         epdf& _epdf() {return *ep;}
 };
 
-/*! \brief Unconditional mpdf, allows using epdf in the role of mpdf. 
+//!DataLink is a connection between epdf and its superordinate (Up)
+//! It is assumed that my val is fully present in "Up"
+class datalink {
+protected:
+        //! Remember how long val should be
+        int valsize;
+        //! Remember how long val of "Up" should be
+        int valupsize;
+        //! val-to-val link, indeces of the upper val
+        ivec v2v_up;
+        public:
+        //! Constructor
+        datalink ( const RV &rv, const RV &rv_up ) :
+                valsize(rv.count()), valupsize(rv_up.count()), v2v_up ( rv.dataind ( rv_up ) )  {it_assert_debug(v2v_up.length()==valsize,"rv is not fully in rv_up");}
+        //! Get val for myself from val of "Up"
+        vec get_val ( const vec &val_up ) {it_assert_debug(valupsize==val_up.length(),"Wrong val_up"); return get_vec ( val_up,v2v_up );}
+        //! Fill val of "Up" by my pieces
+        void fill_val ( vec &val_up, const vec &val ) {it_assert_debug(valsize==val.length(),"Wrong val");it_assert_debug(valupsize==val_up.length(),"Wrong val_up");set_subvector ( val_up, v2v_up, val );}
+};
+
+//!DataLink is a connection between mpdf and its superordinate (Up)
+//! This class links 
+class datalink_mpdf: public datalink {
+protected:
+        //!upper_val-to-local_cond link, indeces of the upper val
+        ivec v2c_up;
+        //!upper_val-to-local_cond link, ideces of the local cond
+        ivec v2c_lo;
+        //!cond-to-cond link, indeces of the upper cond
+        ivec c2c_up;
+        //!cond-to-cond link, indeces of the local cond
+        ivec c2c_lo;
+        //! size of cond
+        int csize;
+public:
+        //! Constructor
+        datalink_mpdf ( const mpdf &Me, const mpdf &Up ) :
+                        datalink ( Me._rv(),Up._rv() ) {
+                //establish v2c connection
+                Me._rvc().dataind ( Up._rv(), v2c_lo, v2c_up );
+                //establish c2c connection
+                Me._rvc().dataind ( Up._rvc(), c2c_lo, c2c_up );
+                csize = Me._rvc().count();
+        }
+        //! Get cond for myself from val and cond of "Up"
+        vec get_cond ( const vec &val_up, const vec &cond_up ) {
+                vec tmp(csize);
+                set_subvector(tmp,v2c_lo,val_up(v2c_up));
+                set_subvector(tmp,c2c_lo,cond_up(c2c_up));
+                return tmp;
+        }
+        //! Fill 
+
+};
+
+/*! \brief Unconditional mpdf, allows using epdf in the role of mpdf.
 
 WARNING: the class does not check validity of the \c ep pointer nor its existence.
@@ -251 +309 @@
 public:
         //!Default constructor
-        mepdf (epdf &em ) :mpdf ( em._rv(),RV() ) {ep=&em;};
-        void condition(const vec &cond){}
-};
-
-//!\brief Abstract composition of pdfs, a base for specific classes 
-class compositepdf{
-        protected:
-                //!Number of mpdfs in the composite
-                int n;
-                //! Elements of composition
-                Array<mpdf*> mpdfs;
-                //! Indeces of rvs in common rv
-                Array<ivec> rvsinrv;
-                //! Indeces of rvc in common rv
-                Array<ivec> irvcs_rv;
-                //! Indeces of common rv in rvc
-                Array<ivec> irv_rvcs;
-        public:
-                compositepdf(Array<mpdf*> A0): n(A0.length()), mpdfs(A0), rvsinrv(n), irvcs_rv(n),irv_rvcs(n){};
-                //! find common rv, flag \param checkoverlap modifies whether overlaps are acceptable
-                RV getrv(bool checkoverlap=false);
-                //! common rvc of all mpdfs is written to rvc
-                void setrvc(const RV &rv, RV &rvc);
-                //! fill all rv*inrv* according to 
-                void setindices(const RV &rv);
+        mepdf ( epdf &em ) :mpdf ( em._rv(),RV() ) {ep=&em;};
+        void condition ( const vec &cond ) {}
+};
+
+//!\brief Abstract composition of pdfs, a base for specific classes
+class compositepdf {
+protected:
+        //!Number of mpdfs in the composite
+        int n;
+        //! Elements of composition
+        Array<mpdf*> mpdfs;
+        //! Indeces of rvs in common rv
+        Array<ivec> rvsinrv;
+        //! Indeces of rvc in common rv
+        Array<ivec> irvcs_rv;
+        //! Indeces of common rv in rvc
+        Array<ivec> irv_rvcs;
+public:
+        compositepdf ( Array<mpdf*> A0 ) : n ( A0.length() ), mpdfs ( A0 ), rvsinrv ( n ), irvcs_rv ( n ),irv_rvcs ( n ) {};
+        //! find common rv, flag \param checkoverlap modifies whether overlaps are acceptable
+        RV getrv ( bool checkoverlap=false );
+        //! common rvc of all mpdfs is written to rvc
+        void setrvc ( const RV &rv, RV &rvc );
+        //! fill all rv*inrv* according to
+        void setindices ( const RV &rv );
 };
 
@@ -327 +385 @@
 
         //!Default constructor
-        BM ( const RV &rv0, double ll0=0,bool evalll0=true ) :rv ( rv0 ), ll ( ll0 ),evalll ( evalll0) {//Fixme: test rv
+        BM ( const RV &rv0, double ll0=0,bool evalll0=true ) :rv ( rv0 ), ll ( ll0 ),evalll ( evalll0 ) {//Fixme: test rv
         };
         //!Copy constructor
-        BM (const BM &B) : rv(B.rv), ll(B.ll), evalll(B.evalll) {}
+        BM ( const BM &B ) : rv ( B.rv ), ll ( B.ll ), evalll ( B.evalll ) {}
 
         /*! \brief Incremental Bayes rule
@@ -337 +395 @@
         virtual void bayes ( const vec &dt ) = 0;
         //! Batch Bayes rule (columns of Dt are observations)
-        virtual void bayesB (const mat &Dt );
+        virtual void bayesB ( const mat &Dt );
         //! Returns a pointer to the epdf representing posterior density on parameters. Use with care!
         virtual const epdf& _epdf() const =0;
@@ -343 +401 @@
         //! Evaluates predictive log-likelihood of the given data record
         //! I.e. marginal likelihood of the data with the posterior integrated out.
-        virtual double logpred(const vec &dt)const{it_error("Not implemented");return 0.0;}
+        virtual double logpred ( const vec &dt ) const{it_error ( "Not implemented" );return 0.0;}
         //! Matrix version of logpred
-        vec logpred_m(const mat &dt)const{vec tmp(dt.cols());for(int i=0;i<dt.cols();i++){tmp(i)=logpred(dt.get_col(i));}return tmp;}
-        
+        vec logpred_m ( const mat &dt ) const{vec tmp ( dt.cols() );for ( int i=0;i<dt.cols();i++ ) {tmp ( i ) =logpred ( dt.get_col ( i ) );}return tmp;}
+
         //!Constructs a predictive density (marginal density on data)
-        virtual epdf* predictor(const RV &rv){it_error("Not implemented");return NULL;};
-        
+        virtual epdf* predictor ( const RV &rv ) {it_error ( "Not implemented" );return NULL;};
+
         //! Destructor for future use;
         virtual ~BM() {};
@@ -356 +414 @@
         //!access function
         double _ll() const {return ll;}
-        //!access function 
-        void set_evalll(bool evl0){evalll=evl0;}
-        
+        //!access function
+        void set_evalll ( bool evl0 ) {evalll=evl0;}
+
         //! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually!
         //! Prototype: BM* _copy_(){BM Tmp*=new Tmp(this*);  return Tmp; }
-        virtual BM* _copy_(bool changerv=false){it_error("function _copy_ not implemented for this BM"); return NULL;};
+        virtual BM* _copy_ ( bool changerv=false ) {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
 };
 

matlab/merger_iter_debug.m

@@ -3 +3 @@
 ndat = 1000; % check if true!
 
+XL = [-1 4];
+YL= XL;
+
 figure(1);
-for it=0:6
+for it=0:20
     figure(it+1)
     subplot(2,2,1);
-    hold off
     si = num2str(it);
-    eval(['plot3(Smp' si '(1,:),Smp' si '(2,:),exp(Mpdf' si '),''r.'')']);
-    if i==0, title('Proposal & mix'); end
-    set(gca,'XLim',[-10,10]);
+    eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),exp(Mpdf' si '))']);
+    title('Proposal density ');
+    set(gca,'XLim',XL);
+    set(gca,'YLim',YL);
         
     subplot(2,2,2);
-    hold off
     eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),exp(lW' si '(1,:)))']);
-%    hold on
-%    eval(['plot3(Smp' si '(1,:),Smp' si '(2,:),w' si ',''r.'')']);
-    if i==0, title('Weighed Smp'); end
+    title('First source'); 
+    set(gca,'XLim',XL);
+    set(gca,'YLim',YL);
     
     subplot(2,2,3);
-    hold off
     eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),exp(lW' si '(2,:)))']);
-%    eval(['plot3(Smp' si '(1,:),Smp' si '(2,:),w' si ',''.'')']);
-%    hold on
-%    eval(['plot3(Smp' si '(1,:),Smp' si '(2,:),w_is_' si ',''r.'')']);
-    if i==0, title('Weighed Smp'); end
+    title('Second source');
+    set(gca,'XLim',XL);
+    set(gca,'YLim',YL);
     
     
     subplot(2,2,4);
     eval(['contour_2(Smp' si '(1,:),Smp' si '(2,:),w' si ')']);
-    if i==0, title('Resampled hist'); end
+    title('Merged density'); 
+    set(gca,'XLim',XL);
+    set(gca,'YLim',YL);
     
-  %  keyboard
 end
 
@@ -43 +44 @@
 figure(it+2);
 M1 = reshape(exp(Res1),Npoints,Npoints);
-contour(XG,YG,M1);
+contour(XG,YG,M1,7);
 
 figure(it+3);
-M2 = reshape(exp(Res2),Npoints,Npoints);
-contour(XG,YG,M2);
-
-figure(it+4);
-L1 = reshape(exp(lf1),Npoints,Npoints);
-contour(XG,YG,L1);
+hold off
+mm = max(max(Res2));
+for i=1:size(Res2,1)
+    M2 = reshape(Res2(i,:),Npoints,Npoints);
+    contour(XG,YG,M2,[0:mm/7:mm]);
+    hold on
+end

tests/CMakeLists.txt

@@ -26 +26 @@
 # BASIC TESTS
 TEST(rv_test)
+TEST(datalink_test)
 TEST(loggers_test)
 

tests/merger_iter_test.cpp

@@ -11 +11 @@
 int main() {
 
+        RNG_randomize();
+        
         RV x ( "{x }","1" );
         RV y ( "{y }","1" );
@@ -16 +18 @@
         RV xy=x; xy.add(y);
         
-//      enorm<fsqmat> f1 ( xy );
-//      enorm<fsqmat> f2 ( x );
-
-        enorm<ldmat> f(xy);
-        f.set_parameters ( "1 1",mat ( "1.2 1; 1 1.2" ) );
-        
-        mpdf* f1=f.condition(x);
-        mpdf* f2=f.condition(y);
-        cout << *(mlnorm<ldmat>*)f1 <<endl;
-        cout << *(mlnorm<ldmat>*)f2 <<endl;
-        
-//      f1.set_parameters ( "1 1",mat ( "1.2 1; 1 1.2" ) );
-//      f2.set_parameters ( "2",mat ( "20" ) );
+        enorm<fsqmat> f1 ( xy );
+        enorm<fsqmat> f2 ( x );
+                
+        f1.set_parameters ( "0 0",mat ( "0.5 0; 0 0.5" ) );
+        f2.set_parameters ( "3",mat ( "0.5" ) );
         
         Array<mpdf* > A ( 2 );
-/*      mepdf A1(f1);
-        mepdf A2(f2);*/
-        A ( 0 ) =f1;
-        A ( 1 ) =f2;
+        mepdf A1(f1);
+        mepdf A2(f2);
+        A ( 0 ) =&A1;
+        A ( 1 ) =&A2;
         
         int Npoints=100;
@@ -48 +42 @@
         merger M ( A );
         enorm<fsqmat> g0(xy);
-        g0.set_parameters(vec("1 1"),mat("100 0; 0 100"));
+        g0.set_parameters(vec("1 1"),mat("1 0; 0 1"));
         
-        M.set_parameters(1.2,200,1);
+        M.set_parameters(1.2,200,2);
         M.merge(&g0);
         
@@ -57 +51 @@
         
         vec Res1 = M.evalpdflog_m(Grid);
-        vec Res2 = MP->evalpdflog_m(Grid);
+        mat Res2 = ((emix*)MP)->evalpdflog_M(Grid);
         
         it_file it("merger_iter_test.it");