Changeset 395

Timestamp:

06/22/09 13:17:49 (15 years ago)

Author:

smidl

Message:

merging works for merger_mx

Files:

• applications/bdmtoolbox/mex/merger_mx.cpp (modified) (2 diffs)
• library/bdm/base/bdmbase.h (modified) (5 diffs)
• library/bdm/stat/emix.h (modified) (3 diffs)
• library/bdm/stat/exp_family.cpp (modified) (3 diffs)
• library/bdm/stat/exp_family.h (modified) (7 diffs)
• library/bdm/stat/merger.cpp (modified) (2 diffs)
• library/bdm/stat/merger.h (modified) (7 diffs)

applications/bdmtoolbox/mex/merger_mx.cpp

@@ -83 +83 @@
         Merger->set_sources(Sources,true); // takes care of deletion of sources
         Merger->merge();
-        
+
+        mxArray* tmp ;
         // Save results
         if (n_output>0){
-                mxArray* tmp = mxCreateStructMatrix(1,1,0,NULL);
+                tmp = mxCreateStructMatrix(1,1,0,NULL);
                 //support
                 Array<vec> &samples=Merger->_Smp()._samples();
@@ -100 +101 @@
                 vec2mxArray(w,fldw);
                 mxReplaceFieldNM(tmp, "weights", fldw);
-                
+
+                // sources
+                                char srcstr[20];
+                for (int i=0;i<Sources.length();i++){
+                        sprintf(srcstr,"source%d",i+1);
+                        vec sll=exp(Sources(i)->evallogcond_m(Merger->_Smp()._samples(),vec(0)));
+
+                        mxArray* fldw=mxCreateDoubleMatrix(1, sll.length(), mxREAL);
+                        vec2mxArray(sll/sum(sll),fldw);
+                        mxReplaceFieldNM(tmp, srcstr, fldw);
+                }               
+
                 output[0] = tmp;
         }

library/bdm/base/bdmbase.h

@@ -268 +268 @@
   //! Compute log-probability of multiple values argument \c val
   virtual vec evallog_m(const mat &Val) const {
-    vec x(Val.cols());
-    for (int i = 0; i < Val.cols(); i++) {x(i) = evallog(Val.get_col(i)) ;}
-    return x;
+          vec x(Val.cols());
+          for (int i = 0; i < Val.cols(); i++) {x(i) = evallog(Val.get_col(i)) ;}
+          return x;
+  }
+  //! Compute log-probability of multiple values argument \c val
+  virtual vec evallog_m(const Array<vec> &Avec) const {
+          vec x(Avec.size());
+          for (int i = 0; i < Avec.size(); i++) {x(i) = evallog(Avec(i)) ;}
+          return x;
   }
   //! Return conditional density on the given RV, the remaining rvs will be in conditioning
@@ -311 +317 @@
   //! Size of the random variable
   int dimension() const {return dim;}
-  //!@}
+  //! Load from structure with elements:
+  //!  \code
+  //! { rv = {class="RV", names=(...),}; // RV describing meaning of random variable
+  //!   // elements of offsprings
+  //! }
+  //! \endcode
+  //!@}
+  void from_setting(const Setting &set){
+          if (set.exists("rv")){
+                  RV* r = UI::build<RV>(set,"rv");
+                  set_rv(*r); 
+                  delete r;
+          }
+  }
 
 };
@@ -368 +387 @@
   //! Matrix version of evallogcond
   virtual vec evallogcond_m(const mat &Dt, const vec &cond) {this->condition(cond); return ep->evallog_m(Dt);};
+  //! Array<vec> version of evallogcond
+  virtual vec evallogcond_m(const Array<vec> &Dt, const vec &cond) {this->condition(cond); return ep->evallog_m(Dt);};
 
   //! \name Access to attributes
@@ -378 +399 @@
   epdf& _epdf() {return *ep;}
   epdf* _e() {return ep;}
+  //! Load from structure with elements:
+  //!  \code
+  //! { rv = {class="RV", names=(...),}; // RV describing meaning of random variable
+  //!   rvc= {class="RV", names=(...),}; // RV describing meaning of random variable in condition
+  //!   // elements of offsprings
+  //! }
+  //! \endcode
+  //!@}
+  void from_setting(const Setting &set){
+          if (set.exists("rv")){
+                  RV* r = UI::build<RV>(set,"rv");
+                  set_rv(*r); 
+                  delete r;
+          }
+          if (set.exists("rvc")){
+                  RV* r = UI::build<RV>(set,"rvc");
+                  set_rvc(*r); 
+                  delete r;
+          }
+  }
   //!@}
 
@@ -570 +611 @@
 public:
         //!Default constructor
+        mepdf(){};
         mepdf ( epdf* em, bool owning_ep0=false ) :mpdf ( ) {ep= em ;owning_ep=owning_ep0;dimc=0;};
         mepdf (const epdf* em ) :mpdf ( ) {ep=const_cast<epdf*>( em );};
         void condition ( const vec &cond ) {}
         ~mepdf(){if (owning_ep) delete ep;}
-};
+  //! Load from structure with elements:
+  //!  \code
+  //! { class = "mepdf",         
+  //!   epdfs = {class="epdfs",...}
+  //! }
+  //! \endcode
+  //!@}
+        void from_setting(const Setting &set){
+                epdf* e = UI::build<epdf>(set,"epdf");
+                ep=     e; 
+                owning_ep=true;
+        }
+};
+UIREGISTER(mepdf);
 
 //!\brief Chain rule of pdfs - abstract part common for mprod and merger. 

library/bdm/stat/emix.h

@@ -235 +235 @@
         mprod (){};
         mprod (Array<mpdf*> mFacs ){set_elements( mFacs );};
-        void set_elements(Array<mpdf*> mFacs ) {
+        void set_elements(Array<mpdf*> mFacs , bool own=false) {
                 
-                set_elements(mFacs);
-                
+                compositepdf::set_elements(mFacs,own);
+                dls.set_size(mFacs.length());
+                epdfs.set_size(mFacs.length());
+                                
                 ep=&dummy;
                 RV rv=getrv ( true );
-                set_rv ( rv );dummy.set_parameters ( rv._dsize() );
+                set_rv ( rv );
+                dummy.set_parameters ( rv._dsize() );
                 setrvc ( ep->_rv(),rvc );
                 // rv and rvc established = > we can link them with mpdfs
@@ -270 +273 @@
                 return res;
         }
+        vec evallogcond_m(const mat &Dt, const vec &cond) {
+                vec tmp(Dt.cols());
+                for(int i=0;i<Dt.cols(); i++){
+                        tmp(i) = evallogcond(Dt.get_col(i),cond);
+                }
+                return tmp;
+        };
+        vec evallogcond_m(const Array<vec> &Dt, const vec &cond) {
+                vec tmp(Dt.length());
+                for(int i=0;i<Dt.length(); i++){
+                        tmp(i) = evallogcond(Dt(i),cond);
+                }
+                return tmp;             
+        };
+
+
         //TODO smarter...
         vec samplecond ( const vec &cond ) {
@@ -294 +313 @@
 
         ~mprod() {};
-};
+        //! Load from structure with elements:
+        //!  \code
+        //! { class='mprod';
+        //!   mpdfs = (..., ...);     // list of mpdfs in the order of chain rule
+        //! }
+        //! \endcode
+        //!@}
+        void from_setting(const Setting &set){
+                Array<mpdf*> Atmp; //temporary Array
+                UI::get(Atmp,set, "mpdfs");
+                set_elements(Atmp,true);
+        }
+        
+};
+UIREGISTER(mprod);
 
 //! Product of independent epdfs. For dependent pdfs, use mprod.

library/bdm/stat/exp_family.cpp

@@ -24 +24 @@
         if ( dimx==1 ) { //same as the following, just quicker.
                 double r = val ( vend ); //last entry!
+                if (r<0) return -1e+100;
                 vec Psi ( nPsi+dimx );
                 Psi ( 0 ) = -1.0;
@@ -34 +35 @@
                 mat Th= reshape ( val ( 0,nPsi*dimx-1 ),nPsi,dimx );
                 fsqmat R ( reshape ( val ( nPsi*dimx,vend ),dimx,dimx ) );
+                double ldetR=R.logdet();
+                if (ldetR) return -1e+100;
                 mat Tmp=concat_vertical ( -eye ( dimx ),Th );
                 fsqmat iR ( dimx );
                 R.inv ( iR );
 
-                return -0.5* ( nu*R.logdet() + trace ( iR.to_mat() *Tmp.T() *V.to_mat() *Tmp ) );
+                return -0.5* ( nu*ldetR + trace ( iR.to_mat() *Tmp.T() *V.to_mat() *Tmp ) );
         }
 }
@@ -187 +190 @@
         int i;
 
+        if (any(val<=0.)) return -1e100;
+        if (any(beta<=0.)) return -1e100;
         for ( i=0; i<dim; i++ ) {
                 res += ( alpha ( i ) - 1 ) *std::log ( val ( i ) ) - beta ( i ) *val ( i );

library/bdm/stat/exp_family.h

@@ -48 +48 @@
                         virtual double evallog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
                         //!Evaluate normalized log-probability
-                        virtual double evallog ( const vec &val ) const {double tmp;tmp= evallog_nn ( val )-lognc();it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); return tmp;}
+                        virtual double evallog ( const vec &val ) const {
+                                double tmp;
+                                tmp= evallog_nn ( val )-lognc();
+                                it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); 
+                                return tmp;}
                         //!Evaluate normalized log-probability for many samples
                         virtual vec evallog ( const mat &Val ) const
@@ -125 +129 @@
                         void set_parameters ( const vec &mu,const sq_T &R );
                         void from_setting(const Setting &root);
+                        void validate() {
+                                it_assert(mu.length()==R.rows(),"parameters mismatch");
+                                dim = mu.length();
+                        }
                         //!@}
 
@@ -225 +233 @@
                         double& _nu()  {return nu;}
                         const double& _nu() const {return nu;}
-                        //!@}
-        };
+                        void from_setting(const Setting &set){
+                                set.lookupValue("nu",nu);
+                                set.lookupValue("dimx",dimx);
+                                mat V;
+                                UI::get(V,set,"V");
+                                set_parameters(dimx, V, nu);
+                                RV* rv=UI::build<RV>(set,"rv");
+                                set_rv(*rv);
+                                delete rv;
+                        }
+                        //!@}
+        };
+        UIREGISTER(egiw);
 
         /*! \brief Dirichlet posterior density
@@ -370 +389 @@
                         vec mean() const {return elem_div ( alpha,beta );}
                         vec variance() const {return elem_div ( alpha,elem_mult ( beta,beta ) ); }
-        };
-
+                        
+                        //! Load from structure with elements:
+                        //!  \code
+                        //! { alpha = [...];         // vector of alpha
+                        //!   beta = [...];          // vector of beta
+                        //!   rv = {class="RV",...}  // description
+                        //! }
+                        //! \endcode
+                        //!@}
+                        void from_setting(const Setting &set){
+                                epdf::from_setting(set); // reads rv
+                                UI::get(alpha,set,"alpha");
+                                UI::get(beta,set,"beta");
+                                validate();
+                        }
+                        void validate(){
+                                it_assert(alpha.length() ==beta.length(), "parameters do not match");
+                                dim =alpha.length();
+                        }
+        };
+UIREGISTER(egamma);
         /*!
          \brief Inverse-Gamma posterior density
@@ -462 +500 @@
                         vec mean() const {return ( high-low ) /2.0;}
                         vec variance() const {return ( pow ( high,2 ) +pow ( low,2 ) +elem_mult ( high,low ) ) /3.0;}
+                        //! Load from structure with elements:
+                        //!  \code
+                        //! { high = [...];          // vector of upper bounds
+                        //!   low = [...];           // vector of lower bounds
+                        //!   rv = {class="RV",...}  // description of RV
+                        //! }
+                        //! \endcode
+                        //!@}
+                        void from_setting(const Setting &set){
+                                epdf::from_setting(set); // reads rv and rvc
+                                UI::get(high,set,"high");
+                                UI::get(low,set,"low");
+                        }
         };
 
@@ -650 +701 @@
                         void set_parameters ( double k, const vec &beta0 );
                         void condition ( const vec &val ) {_beta=k/val;};
-        };
-
+                        //! Load from structure with elements:
+                        //!  \code
+                        //! { alpha = [...];         // vector of alpha
+                        //!   k = 1.1;               // multiplicative constant k
+                        //!   rv = {class="RV",...}  // description of RV
+                        //!   rvc = {class="RV",...} // description of RV in condition
+                        //! }
+                        //! \endcode
+                        //!@}
+                        void from_setting(const Setting &set){
+                                mpdf::from_setting(set); // reads rv and rvc
+                                vec betatmp; // ugly but necessary
+                                UI::get(betatmp,set,"beta");
+                                set.lookupValue("k",k);
+                                set_parameters(k,betatmp);
+                        }
+        };
+        UIREGISTER(mgamma);
+        
         /*!
         \brief  Inverse-Gamma random walk
@@ -1088 +1156 @@
                 mu = mu0;
                 R = R0;
-                dim = mu0.length();
-        };
-
-        template<class sq_T>
-        void enorm<sq_T>::from_setting(const Setting &root){
-                vec mu;
-                UI::get(mu,root,"mu");
-                mat R;
-                UI::get(R,root,"R");
-                set_parameters(mu,R);
+                validate();
+        };
+
+        template<class sq_T>
+        void enorm<sq_T>::from_setting(const Setting &set){
+                epdf::from_setting(set); //reads rv
                 
-                RV* r = UI::build<RV>(root,"rv");
-                set_rv(*r); 
-                delete r;
+                UI::get(mu,set,"mu");
+                mat Rtmp;// necessary for conversion
+                UI::get(Rtmp,set,"R");
+                R=Rtmp; // conversion
+                validate();
         }
 

library/bdm/stat/merger.cpp

@@ -134 +134 @@
                                 {
                                         // no need for conditioning or marginalization
-                                        for ( int j=0;j<Npoints; j++ )   // Smp is Array<> => for cycle
-                                        {
-                                                lw_src ( j ) =mpdfs ( i )->_epdf().evallog ( Smp ( j ) );
-                                        }
+                                        lw_src = mpdfs ( i )->_epdf().evallog_m ( Smp  );
                                 }
                                 else
@@ -219 +216 @@
                         // ==== stopping rule ===
                         niter++;
-                        converged = ( niter>40 );
+                        converged = ( niter>stop_niter );
                 }
                 delete Mpred;

library/bdm/stat/merger.h

@@ -132 +132 @@
                         }
                         // fill samples
-                        int act_dim=0; //active dimension
                         for (int i = 0; i < Npoints; i++) {
                                 // copy 
@@ -139 +138 @@
                                 for (int j = 0;j < dim;j++) {
                                         if (ind (j) == gridsize (j) - 1) { //j-th index is full
-//                                              ind (j) = 0; //shift back
+                                                ind (j) = 0; //shift back
                                                 smpi(j) = XYZ(j)(0);
                                                 
-//                                              ind (j + 1) ++; //increase the next dimension;
-                                                smpi(j+1) += steps(j+1);
+                                                if (i<Npoints-1) {
+                                                        ind (j + 1) ++; //increase the next dimension;
+                                                        smpi(j+1) += steps(j+1);
+                                                        break;
+                                                }
                                                 
-                                                if (ind (j + 1) < gridsize (j + 1) - 1) break;
                                         } else {
-//                                              ind (j) ++; 
+                                                ind (j) ++; 
                                                 smpi(j) +=steps(j);
                                                 break;
@@ -185 +186 @@
 
                 //!Merge given sources in given points
-                void merge () {
+                virtual void merge () {
                         validate();
 
@@ -205 +206 @@
                                 }
 
-                                vec wtmp = exp (merge_points (lW));
+                                vec w_nn=merge_points (lW);
+                                vec wtmp = exp (w_nn-max(w_nn));
                                 //renormalize
                                 eSmp._w() = wtmp / sum (wtmp);
@@ -277 +279 @@
                                 }
                         }
+                        if (set.exists("dbg_file")){ 
+                                string dbg_file;
+                                UI::get<string> (dbg_file, set, "dbg_file");
+                                set_debug_file(dbg_file);
+                        }
+                        //validate() - not used
                 }
 
@@ -295 +303 @@
                 //!Number of components in a mixture
                 int Ncoms;
-                //! coefficient of resampling
+                //! coefficient of resampling [0,1]
                 double effss_coef;
+                //! stop after niter iterations
+                int stop_niter;
 
         public:
@@ -338 +348 @@
                 emix* proposal() {emix* tmp = Mix.epredictor(); tmp->set_rv (rv); return tmp;}
                 //! from_settings
-                void from_settings(const Setting& set){
+                void from_setting(const Setting& set){
                         merger_base::from_setting(set);
                         set.lookupValue("ncoms",Ncoms);
+                        set.lookupValue("effss_coef",effss_coef);
+                        set.lookupValue("stop_niter",stop_niter);
                 }