Index: library/bdm/estim/particles.cpp
===================================================================
--- library/bdm/estim/particles.cpp (revision 487)
+++ library/bdm/estim/particles.cpp (revision 638)
@@ -5,33 +5,48 @@
 using std::endl;
 
-void PF::bayes ( const vec &dt ) {
-	int i;
-	vec lls ( n );
-	ivec ind;
-	double mlls = -std::numeric_limits<double>::infinity(), sum = 0.0;
+void PF::bayes_gensmp(){
+	#pragma parallel for
+	for (int i = 0; i < n; i++ ) {
+		_samples ( i ) = par->samplecond ( _samples ( i ) );
+		lls(i) = 0;
+	}
+}
 
-	for ( i = 0; i < n; i++ ) {
-		//generate new samples from paramater evolution model;
-		vec old_smp=_samples ( i );
-		_samples ( i ) = par->samplecond ( old_smp );
-		lls ( i ) = par->evallogcond ( _samples ( i ), old_smp );
-		lls ( i ) *= obs->evallogcond ( dt, _samples ( i ) );
-
-		if ( lls ( i ) > mlls ) mlls = lls ( i ); //find maximum
-	}
-
+void PF::bayes_weights(){
+// 
+	double mlls=max(lls);
 	// compute weights
-	for ( i = 0; i < n; i++ ) {
+	for (int  i = 0; i < n; i++ ) {
 		_w ( i ) *= exp ( lls ( i ) - mlls ); // multiply w by likelihood
 	}
 
 	//renormalize
+	double sw=sum(_w);
+	if (!isfinite(sw)) {
+		for (int i=0;i<n;i++){
+			if (!isfinite(_w(i))) {_w(i)=0;}
+		}
+		sw = sum(_w);
+		if (!isfinite(sw)) {
+			bdm_error("Particle filter is lost; no particle is good enough.");
+		}
+	}
+	_w /= sw;
+}
+
+void PF::bayes ( const vec &dt ) {
+	int i;
+	// generate samples - time step
+	bayes_gensmp();
+	// weight them - data step
 	for ( i = 0; i < n; i++ ) {
-		sum += _w ( i );
-	};
+		lls ( i ) += obs->evallogcond ( dt, _samples ( i ) ); //+= because lls may have something from gensmp!
+	}
 
-	_w ( i ) /= sum; //?
+	bayes_weights();
 
-	ind = est.resample ( resmethod );
+	if (do_resampling()) {
+		est.resample ( resmethod);
+	}
 
 }
@@ -45,4 +60,37 @@
 // }
 
+void MPF::bayes ( const vec &dt ) {	
+	// follows PF::bayes in most places!!!
+	
+	int i;
+	int n=pf->__w().length();
+	vec &lls = pf->_lls();
+	
+	// generate samples - time step
+	pf->bayes_gensmp();
+	// weight them - data step
+	#pragma parallel for
+	for ( i = 0; i < n; i++ ) {
+		BMs(i) -> condition(pf->posterior()._sample(i));
+		BMs(i) -> bayes(dt);
+		lls ( i ) += BMs(i)->_ll();
+	}
+
+	pf->bayes_weights();
+	
+	ivec ind;
+	if (pf->do_resampling()) {
+		pf->resample ( ind);
+		
+		#pragma omp parallel for
+		for ( i = 0; i < n; i++ ) {
+			if ( ind ( i ) != i ) {//replace the current Bm by a new one
+				delete BMs ( i );
+				BMs ( i ) = BMs ( ind ( i ) )->_copy_(); //copy constructor
+			}
+		};
+	}
+};
+
 
 }
Index: library/bdm/estim/particles.h
===================================================================
--- library/bdm/estim/particles.h (revision 565)
+++ library/bdm/estim/particles.h (revision 638)
@@ -36,11 +36,16 @@
 	Array<vec> &_samples;
 	//! Parameter evolution model
-	mpdf *par;
+	shared_ptr<mpdf> par;
 	//! Observation model
-	mpdf *obs;
-
+	shared_ptr<mpdf> obs;
+	//! internal structure storing loglikelihood of predictions
+	vec lls;
+	
 	//! which resampling method will be used
 	RESAMPLING_METHOD resmethod;
-
+	//! resampling threshold; in this case its meaning is minimum ratio of active particles
+	//! For example, for 0.5 resampling is performed when the numebr of active aprticles drops belo 50%.
+	double res_threshold;
+	
 	//! \name Options
 	//!@{
@@ -58,19 +63,28 @@
 		LIDs.set_size ( 5 );
 	};
-	/*	PF ( mpdf *par0, mpdf *obs0, epdf *epdf0, int n0 ) :
-				est ( ),_w ( est._w() ),_samples ( est._samples() ),opt_L_smp(false), opt_L_wei(false)
-		{ set_parameters ( par0,obs0,n0 ); set_statistics ( ones ( n0 ),epdf0 ); };*/
-	void set_parameters ( mpdf *par0, mpdf *obs0, int n0, RESAMPLING_METHOD rm = SYSTEMATIC ) {
+	
+	void set_parameters (int n0, double res_th0=0.5, RESAMPLING_METHOD rm = SYSTEMATIC ) {
+		n = n0;
+		res_threshold = res_th0;
+		resmethod = rm;
+	};
+	void set_model ( shared_ptr<mpdf> par0, shared_ptr<mpdf> obs0) {
 		par = par0;
 		obs = obs0;
-		n = n0;
-		resmethod = rm;
+		// set values for posterior
+		est.set_rv(par->_rv());
 	};
 	void set_statistics ( const vec w0, const epdf &epdf0 ) {
 		est.set_statistics ( w0, epdf0 );
 	};
+	void set_statistics ( const eEmp &epdf0 ) {
+		bdm_assert_debug(epdf0._rv().equal(par->_rv()),"Incompatibel input");
+		est=epdf0;
+	};
 	//!@}
 	//! Set posterior density by sampling from epdf0
-//	void set_est ( const epdf &epdf0 );
+	//! Extends original BM::set_options by two more options:
+	//! \li logweights - meaning that all weightes will be logged
+	//! \li logsamples - all samples will be also logged
 	void set_options ( const string &opt ) {
 		BM::set_options ( opt );
@@ -78,76 +92,163 @@
 		opt_L_smp = ( opt.find ( "logsamples" ) != string::npos );
 	}
+	//! bayes I - generate samples and add their weights to lls
+	virtual void bayes_gensmp();
+	//! bayes II - compute weights of the 
+	virtual void bayes_weights();
+	//! important part of particle filtering - decide if it is time to perform resampling
+	virtual bool do_resampling(){ 	
+		double eff = 1.0 / ( _w * _w );
+		return eff < ( res_threshold*n );
+	}
 	void bayes ( const vec &dt );
 	//!access function
-	vec* __w() {
-		return &_w;
+	vec& __w() { return _w; }
+	//!access function
+	vec& _lls() { return lls; }
+	RESAMPLING_METHOD _resmethod() const { return resmethod; }
+	//!access function
+	const eEmp& posterior() const {return est;}
+	
+	/*! configuration structure for basic PF
+	\code
+	parameter_pdf   = mpdf_class;         // parameter evolution pdf
+	observation_pdf = mpdf_class;         // observation pdf
+	prior           = epdf_class;         // prior probability density
+	--- optional ---
+	n               = 10;                 // number of particles
+	resmethod       = 'systematic', or 'multinomial', or 'stratified'
+										  // resampling method
+	res_threshold   = 0.5;                // resample when active particles drop below 50%
+	\endcode
+	*/
+	void from_setting(const Setting &set){
+		par = UI::build<mpdf>(set,"parameter_pdf",UI::compulsory);
+		obs = UI::build<mpdf>(set,"observation_pdf",UI::compulsory);
+		
+		prior_from_set(set);
+		resmethod_from_set(set);
+		// set resampling method
+		//set drv
+		//find potential input - what remains in rvc when we subtract rv
+		RV u = par->_rvc().remove_time().subt( par->_rv() ); 
+		//find potential input - what remains in rvc when we subtract x_t
+		RV obs_u = obs->_rvc().remove_time().subt( par->_rv() ); 
+		
+		u.add(obs_u); // join both u, and check if they do not overlap
+		
+		set_drv(concat(obs->_rv(),u) );
+	}
+	//! auxiliary function reading parameter 'resmethod' from configuration file
+	void resmethod_from_set(const Setting &set){
+		string resmeth;
+		if (UI::get(resmeth,set,"resmethod",UI::optional)){
+			if (resmeth=="systematic") {
+				resmethod= SYSTEMATIC;
+			} else  {
+				if (resmeth=="multinomial"){
+					resmethod=MULTINOMIAL;
+				} else {
+					if (resmeth=="stratified"){
+						resmethod= STRATIFIED;
+					} else {
+						bdm_error("Unknown resampling method");
+					}
+				}
+			}
+		} else {
+			resmethod=SYSTEMATIC;
+		};
+		if(!UI::get(res_threshold, set, "res_threshold", UI::optional)){
+			res_threshold=0.5;
+		}
+	}
+	//! load prior information from set and set internal structures accordingly
+	void prior_from_set(const Setting & set){
+		shared_ptr<epdf> pri = UI::build<epdf>(set,"prior",UI::compulsory);
+		
+		eEmp *test_emp=dynamic_cast<eEmp*>(&(*pri));
+		if (test_emp) { // given pdf is sampled
+			est=*test_emp;
+		} else {
+			int n;
+			if (!UI::get(n,set,"n",UI::optional)){n=10;}
+			// sample from prior
+			set_statistics(ones(n)/n, *pri);
+		}
+		//validate();
+	}
+	
+	void validate(){
+		n=_w.length();
+		lls=zeros(n);
+		if (par->_rv()._dsize()>0) {
+			bdm_assert(par->_rv()._dsize()==est.dimension(),"Mismatch of RV and dimension of posterior" );
+		}
+	}
+	//! resample posterior density (from outside - see MPF)
+	void resample(ivec &ind){
+		est.resample(ind,resmethod);
 	}
 };
+UIREGISTER(PF);
 
 /*!
 \brief Marginalized Particle filter
 
-Trivial version: proposal = parameter evolution, observation model is not used. (it is assumed to be part of BM).
+A composition of particle filter with exact (or almost exact) bayesian models (BMs).
+The Bayesian models provide marginalized predictive density. Internaly this is achieved by virtual class MPFmpdf.
 */
 
-template<class BM_T>
-class MPF : public PF {
-	Array<BM_T*> BMs;
+class MPF : public BM  {
+    shared_ptr<PF> pf;
+	Array<BM*> BMs;
 
 	//! internal class for MPDF providing composition of eEmp with external components
 
 	class mpfepdf : public epdf  {
-	protected:
-		eEmp &E;
-		vec &_w;
-		Array<const epdf*> Coms;
+		shared_ptr<PF> &pf;
+		Array<BM*> &BMs;
 	public:
-		mpfepdf ( eEmp &E0 ) :
-				epdf ( ), E ( E0 ),  _w ( E._w() ),
-				Coms ( _w.length() ) {
-		};
-		//! read statistics from MPF
-		void read_statistics ( Array<BM_T*> &A ) {
-			dim = E.dimension() + A ( 0 )->posterior().dimension();
-			for ( int i = 0; i < _w.length() ; i++ ) {
-				Coms ( i ) = &(A ( i )->posterior());
+		mpfepdf (shared_ptr<PF> &pf0, Array<BM*> &BMs0): epdf(), pf(pf0), BMs(BMs0) { };
+		//! a variant of set parameters - this time, parameters are read from BMs and pf
+		void read_parameters(){
+			rv = concat(pf->posterior()._rv(), BMs(0)->posterior()._rv());
+			dim = pf->posterior().dimension() + BMs(0)->posterior().dimension();
+			bdm_assert_debug(dim == rv._dsize(), "Wrong name ");
+		}
+		vec mean() const {
+			const vec &w = pf->posterior()._w();
+			vec pom = zeros ( BMs(0)->posterior ().dimension() );
+			//compute mean of BMs
+			for ( int i = 0; i < w.length(); i++ ) {
+				pom += BMs ( i )->posterior().mean() * w ( i );
 			}
-		}
-		//! needed in resampling
-		void set_elements ( int &i, double wi, const epdf* ep ) {
-			_w ( i ) = wi;
-			Coms ( i ) = ep;
-		};
-
-		void set_parameters ( int n ) {
-			E.set_parameters ( n, false );
-			Coms.set_length ( n );
-		}
-		vec mean() const {
-			// ugly
-			vec pom = zeros ( Coms ( 0 )->dimension() );
-			for ( int i = 0; i < _w.length(); i++ ) {
-				pom += Coms ( i )->mean() * _w ( i );
+			return concat ( pf->posterior().mean(), pom );
+		}
+		vec variance() const {
+			const vec &w = pf->posterior()._w();
+			
+			vec pom = zeros ( BMs(0)->posterior ().dimension() );
+			vec pom2 = zeros ( BMs(0)->posterior ().dimension() );
+			vec mea;
+			
+			for ( int i = 0; i < w.length(); i++ ) {
+				// save current mean 
+				mea = BMs ( i )->posterior().mean();
+				pom += mea * w ( i );
+				//compute variance
+				pom2 += ( BMs ( i )->posterior().variance() + pow ( mea, 2 ) ) * w ( i );
 			}
-			return concat ( E.mean(), pom );
-		}
-		vec variance() const {
-			// ugly
-			vec pom = zeros ( Coms ( 0 )->dimension() );
-			vec pom2 = zeros ( Coms ( 0 )->dimension() );
-			for ( int i = 0; i < _w.length(); i++ ) {
-				pom += Coms ( i )->mean() * _w ( i );
-				pom2 += ( Coms ( i )->variance() + pow ( Coms ( i )->mean(), 2 ) ) * _w ( i );
-			}
-			return concat ( E.variance(), pom2 - pow ( pom, 2 ) );
-		}
+			return concat ( pf->posterior().variance(), pom2 - pow ( pom, 2 ) );
+		}
+		
 		void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const {
 			//bounds on particles
 			vec lbp;
 			vec ubp;
-			E.qbounds ( lbp, ubp );
+			pf->posterior().qbounds ( lbp, ubp );
 
 			//bounds on Components
-			int dimC = Coms ( 0 )->dimension();
+			int dimC = BMs ( 0 )->posterior().dimension();
 			int j;
 			// temporary
@@ -160,7 +261,8 @@
 			Ubc = -std::numeric_limits<double>::infinity();
 
-			for ( int i = 0; i < _w.length(); i++ ) {
+			for ( int i = 0; i < BMs.length(); i++ ) {
 				// check Coms
-				Coms ( i )->qbounds ( lbc, ubc );
+				BMs ( i )->posterior().qbounds ( lbc, ubc );
+				//save either minima or maxima
 				for ( j = 0; j < dimC; j++ ) {
 					if ( lbc ( j ) < Lbc ( j ) ) {
@@ -195,21 +297,20 @@
 public:
 	//! Default constructor.
-	MPF () : PF (), jest ( est ) {};
-	void set_parameters ( mpdf *par0, mpdf *obs0, int n0, RESAMPLING_METHOD rm = SYSTEMATIC ) {
-		PF::set_parameters ( par0, obs0, n0, rm );
-		jest.set_parameters ( n0 );//duplication of rm
+	MPF () :  jest (pf,BMs) {};
+	void set_parameters ( shared_ptr<mpdf> par0, shared_ptr<mpdf> obs0, int n0, RESAMPLING_METHOD rm = SYSTEMATIC ) {
+		pf->set_model ( par0, obs0); 
+		pf->set_parameters(n0, rm );
 		BMs.set_length ( n0 );
 	}
-	void set_statistics ( const epdf &epdf0, const BM_T* BMcond0 ) {
-
-		PF::set_statistics ( ones ( n ) / n, epdf0 );
+	void set_BM ( const BM &BMcond0 ) {
+
+		int n=pf->__w().length();
+		BMs.set_length(n);
 		// copy
+		//BMcond0 .condition ( pf->posterior()._sample ( 0 ) );
 		for ( int i = 0; i < n; i++ ) {
-			BMs ( i ) = new BM_T ( *BMcond0 );
-			BMs ( i )->condition ( _samples ( i ) );
-		}
-
-		jest.read_statistics ( BMs );
-		//options
+			BMs ( i ) = BMcond0._copy_();
+			BMs ( i )->condition ( pf->posterior()._sample ( i ) );
+		}
 	};
 
@@ -218,13 +319,8 @@
 		return jest;
 	}
-	//! Set postrior of \c rvc to samples from epdf0. Statistics of BMs are not re-computed! Use only for initialization!
-	/*	void set_est ( const epdf& epdf0 ) {
-			PF::set_est ( epdf0 );  // sample params in condition
-			// copy conditions to BMs
-
-			for ( int i=0;i<n;i++ ) {BMs(i)->condition ( _samples ( i ) );}
-		}*/
+	//! Extends options understood by BM::set_options by option 
+	//! \li logmeans - meaning 
 	void set_options ( const string &opt ) {
-		PF::set_options ( opt );
+		BM::set_options(opt);
 		opt_L_mea = ( opt.find ( "logmeans" ) != string::npos );
 	}
@@ -234,63 +330,50 @@
 		return BMs ( i );
 	}
+	
+	/*! configuration structure for basic PF
+	\code
+	BM              = BM_class;           // Bayesian filtr for analytical part of the model
+	parameter_pdf   = mpdf_class;         // transitional pdf for non-parametric part of the model
+	prior           = epdf_class;         // prior probability density
+	--- optional ---
+	n               = 10;                 // number of particles
+	resmethod       = 'systematic', or 'multinomial', or 'stratified'
+										  // resampling method
+	\endcode
+	*/	
+	void from_setting(const Setting &set){
+		shared_ptr<mpdf> par = UI::build<mpdf>(set,"parameter_pdf",UI::compulsory);
+		shared_ptr<mpdf> obs= new mpdf(); // not used!!
+
+		pf = new PF;
+		// rpior must be set before BM
+		pf->prior_from_set(set);
+		pf->resmethod_from_set(set);
+		pf->set_model(par,obs);
+		
+		shared_ptr<BM> BM0 =UI::build<BM>(set,"BM",UI::compulsory);
+		set_BM(*BM0);
+		
+		string opt;
+		if (UI::get(opt,set,"options",UI::optional)){
+			set_options(opt);
+		}
+		//set drv
+		//find potential input - what remains in rvc when we subtract rv
+		RV u = par->_rvc().remove_time().subt( par->_rv() ); 		
+		set_drv(concat(BM0->_drv(),u) );
+		validate();
+	}
+	void validate(){
+		try{
+		pf->validate();
+		} catch (std::exception &e){
+			throw UIException("Error in PF part of MPF:");
+		}
+		jest.read_parameters();
+	}
+	
 };
-
-template<class BM_T>
-void MPF<BM_T>::bayes ( const vec &dt ) {
-	int i;
-	vec lls ( n );
-	vec llsP ( n );
-	ivec ind;
-	double mlls = -std::numeric_limits<double>::infinity();
-
-#pragma omp parallel for
-	for ( i = 0; i < n; i++ ) {
-		//generate new samples from paramater evolution model;
-		vec old_smp=_samples ( i );
-		_samples ( i ) = par->samplecond ( old_smp );
-		llsP ( i ) = par->evallogcond ( _samples ( i ), old_smp );
-		BMs ( i )->condition ( _samples ( i ) );
-		BMs ( i )->bayes ( dt );
-		lls ( i ) = BMs ( i )->_ll(); // lls above is also in proposal her must be lls(i) =, not +=!!
-		if ( lls ( i ) > mlls ) mlls = lls ( i ); //find maximum likelihood (for numerical stability)
-	}
-
-	double sum_w = 0.0;
-	// compute weights
-#pragma omp parallel for
-	for ( i = 0; i < n; i++ ) {
-		_w ( i ) *= exp ( lls ( i ) - mlls ); // multiply w by likelihood
-		sum_w += _w ( i );
-	}
-
-	if ( sum_w  > 0.0 ) {
-		_w /= sum_w; //?
-	} else {
-		cout << "sum(w)==0" << endl;
-	}
-
-
-	double eff = 1.0 / ( _w * _w );
-	if ( eff < ( 0.3*n ) ) {
-		ind = est.resample ( resmethod );
-		// Resample Bms!
-
-#pragma omp parallel for
-		for ( i = 0; i < n; i++ ) {
-			if ( ind ( i ) != i ) {//replace the current Bm by a new one
-				//fixme this would require new assignment operator
-				// *Bms[i] = *Bms[ind ( i ) ];
-
-				// poor-man's solution: replicate constructor here
-				// copied from MPF::MPF
-				delete BMs ( i );
-				BMs ( i ) = new BM_T ( *BMs ( ind ( i ) ) ); //copy constructor
-				const epdf& pom = BMs ( i )->posterior();
-				jest.set_elements ( i, 1.0 / n, &pom );
-			}
-		};
-		cout << '.';
-	}
-}
+UIREGISTER(MPF);
 
 }
