Changeset 192 for bdm/stat/emix.h

Timestamp:

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

Author:

smidl

Message:

modification of datalinks and switch mprod and merger to use them

Files:

• bdm/stat/emix.h (modified) (4 diffs)

bdm/stat/emix.h

@@ -35 +35 @@
  */
 class mratio: public mpdf {
-        protected:
+protected:
         //! Nominator in the form of mpdf
-                const epdf* nom;
+        const epdf* nom;
         //!Denominator in the form of epdf
-                epdf* den;
-        //!flag for destructor 
-                bool destroynom;
-        public:
-        //!Default constructor. By default, the given epdf is not copied! 
+        epdf* den;
+        //!flag for destructor
+        bool destroynom;
+public:
+        //!Default constructor. By default, the given epdf is not copied!
         //! It is assumed that this function will be used only temporarily.
-                mratio(const epdf* nom0, const RV &rv, bool copy=false):mpdf(rv,RV()){
-                        if (copy){
+        mratio ( const epdf* nom0, const RV &rv, bool copy=false ) :mpdf ( rv,RV() ) {
+                if ( copy ) {
 //                      nom = nom0->_copy_();
-                                destroynom=true;
-                        } else {
-                                nom = nom0;
-                                destroynom = false;
-                        }
-                        rvc = nom->_rv().subt(rv);
-                        it_assert_debug(rvc.length()>0,"Makes no sense to use this object!");
-                        den = nom->marginal(rvc);
-                };
-                double evalcond(const vec &val, const vec &cond) {
-                        return exp(nom->evalpdflog(concat(val,cond)) - den->evalpdflog(cond));
-                }
+                        destroynom=true;
+                }
+                else {
+                        nom = nom0;
+                        destroynom = false;
+                }
+                rvc = nom->_rv().subt ( rv );
+                it_assert_debug ( rvc.length() >0,"Makes no sense to use this object!" );
+                den = nom->marginal ( rvc );
+        };
+        double evalcond ( const vec &val, const vec &cond ) {
+                return exp ( nom->evalpdflog ( concat ( val,cond ) ) - den->evalpdflog ( cond ) );
+        }
         //! Object takes ownership of nom and will destroy it
-                void ownnom(){destroynom=true;}
+        void ownnom() {destroynom=true;}
         //! Default destructor
-                ~mratio(){delete den; if (destroynom) {delete nom;}}
+        ~mratio() {delete den; if ( destroynom ) {delete nom;}}
 };
 
@@ -100 +101 @@
                 int i;
                 double sum = 0.0;
-                for ( i = 0;i < w.length();i++ ) {sum += w ( i ) * exp(Coms ( i )->evalpdflog ( val ));}
+                for ( i = 0;i < w.length();i++ ) {sum += w ( i ) * exp ( Coms ( i )->evalpdflog ( val ) );}
                 return log ( sum );
         };
         vec evalpdflog_m ( const mat &Val ) const {
-                vec x=ones(Val.cols());
-                vec y(Val.cols());
-                for (int i = 0; i < w.length(); i++ ) {
-                        y = w ( i )*exp(Coms ( i )->evalpdflog_m ( Val ) );
-                        elem_mult_inplace(y,x); //result is in x
-                }
-                return log(x);
+                vec x=ones ( Val.cols() );
+                vec y ( Val.cols() );
+                for ( int i = 0; i < w.length(); i++ ) {
+                        y = w ( i ) *exp ( Coms ( i )->evalpdflog_m ( Val ) );
+                        elem_mult_inplace ( y,x ); //result is in x
+                }
+                return log ( x );
         };
         mat evalpdflog_M ( const mat &Val ) const {
-                mat X(w.length(), Val.cols());
-                for (int i = 0; i < w.length(); i++ ) {
-                        X.set_row(i, w ( i )*exp(Coms ( i )->evalpdflog_m ( Val ) ));
+                mat X ( w.length(), Val.cols() );
+                for ( int i = 0; i < w.length(); i++ ) {
+                        X.set_row ( i, w ( i ) *exp ( Coms ( i )->evalpdflog_m ( Val ) ) );
                 }
                 return X;
@@ -143 +144 @@
         //! pointers to epdfs - shortcut to mpdfs()._epdf()
         Array<epdf*> epdfs;
-        //! Indeces of rvc in common rvc
-        Array<ivec> irvcs_rvc;
-        //! Indeces of common rvc in rvcs
-        Array<ivec> irvc_rvcs;
+        //! Data link for each mpdfs
+        Array<datalink_m2m*> dls;
 public:
         /*!\brief Constructor from list of mFacs,
         */
-        mprod ( Array<mpdf*> mFacs ) : compositepdf ( mFacs ), mpdf ( getrv ( true ),RV() ), epdfs ( n ), irvcs_rvc ( n ), irvc_rvcs ( n ) {
+        mprod ( Array<mpdf*> mFacs ) : compositepdf ( mFacs ), mpdf ( getrv ( true ),RV() ), epdfs ( n ), dls ( n ) {
                 setrvc ( rv,rvc );
-                setindices ( rv );
+                // rv and rvc established = > we can link them with mpdfs
                 for ( int i = 0;i < n;i++ ) {
-                        // establish link between rvc and rvcs
-                        rvc.dataind ( mpdfs ( i )->_rvc(), irvc_rvcs ( i ), irvcs_rvc ( i ) );
+                        dls ( i ) = new datalink_m2m ( mpdfs ( i )->_rv(), mpdfs(i)->_rvc(), rv, rvc );
                 }
 
@@ -166 +164 @@
                 int i;
                 double res = 0.0;
-                vec condi;
                 for ( i = n - 1;i >= 0;i-- ) {
-                        if ( ( irvcs_rvc ( i ).length() > 0 ) || ( irvcs_rv ( i ).length() > 0 ) ) {
-                                condi = zeros ( irvcs_rvc ( i ).length() +irvcs_rv ( i ).length() );
-                                //copy part of the rvc into condi
-                                set_subvector ( condi, irvcs_rvc ( i ), get_vec ( cond,irvc_rvcs ( i ) ) );
-                                //copy part of the rv into condi
-                                set_subvector ( condi, irvcs_rv ( i ), get_vec ( val,irv_rvcs ( i ) ) );
-                                mpdfs ( i )->condition ( condi );
+                        if ( mpdfs(i)->_rvc().count() >0) {
+                                mpdfs ( i )->condition ( dls ( i )->get_cond ( val,cond ) );
                         }
                         // add logarithms
-                        res += epdfs ( i )->evalpdflog ( val ( rvsinrv ( i ) ) );
+                        res += epdfs ( i )->evalpdflog ( dls ( i )->get_val ( val ) );
                 }
                 return exp ( res );
         }
         vec samplecond ( const vec &cond, double &ll ) {
-                vec smp=zeros ( rv.count() );
-                vec condi;
+                //! Ugly hack to help to discover if mpfs are not in proper order. Correct solution = check that explicitely.
+                vec smp= std::numeric_limits<double>::infinity() * ones ( rv.count() );
                 vec smpi;
                 ll = 0;
+                // Hard assumption here!!! We are going backwards, to assure that samples that are needed from smp are already generated!
                 for ( int i = ( n - 1 );i >= 0;i-- ) {
-                        if ( ( irvcs_rvc ( i ).length() > 0 ) || ( irvcs_rv ( i ).length() > 0 ) ) {
-                                condi = zeros ( irvcs_rv ( i ).length() + irvcs_rvc ( i ).length() );
-                                // copy data in condition
-                                set_subvector ( condi,irvcs_rvc ( i ), get_vec ( cond, irvc_rvcs ( i ) ) );
-                                // copy data in already generated sample
-                                set_subvector ( condi,irvcs_rv ( i ), get_vec ( smp,irv_rvcs ( i ) ) );
-
-                                mpdfs ( i )->condition ( condi );
+                        if ( mpdfs(i)->_rvc().count() ) {
+                                mpdfs ( i )->condition ( dls ( i )->get_cond ( smp ,cond ) ); // smp is val here!!
                         }
                         smpi = epdfs ( i )->sample();
                         // copy contribution of this pdf into smp
-                        set_subvector ( smp,rvsinrv ( i ), smpi );
+                        dls(i)->fill_val(smp, smpi);
                         // add ith likelihood contribution
                         ll+=epdfs ( i )->evalpdflog ( smpi );