Index: library/bdm/estim/particles.h
===================================================================
--- library/bdm/estim/particles.h (revision 487)
+++ library/bdm/estim/particles.h (revision 488)
@@ -67,5 +67,5 @@
 		resmethod = rm;
 	};
-	void set_statistics ( const vec w0, epdf *epdf0 ) {
+	void set_statistics ( const vec w0, const epdf &epdf0 ) {
 		est.set_statistics ( w0, epdf0 );
 	};
@@ -201,5 +201,5 @@
 		BMs.set_length ( n0 );
 	}
-	void set_statistics ( epdf *epdf0, const BM_T* BMcond0 ) {
+	void set_statistics ( const epdf &epdf0, const BM_T* BMcond0 ) {
 
 		PF::set_statistics ( ones ( n ) / n, epdf0 );
Index: library/bdm/base/bdmbase.h
===================================================================
--- library/bdm/base/bdmbase.h (revision 487)
+++ library/bdm/base/bdmbase.h (revision 488)
@@ -25,5 +25,6 @@
 using namespace std;
 
-namespace bdm {
+namespace bdm
+{
 
 typedef std::map<string, int> RVmap;
@@ -32,14 +33,15 @@
 
 //! Structure of RV, i.e. RVs expanded into a flat list of IDs, used for debugging.
-class str {
-public:
-	//! vector id ids (non-unique!)
-	ivec ids;
-	//! vector of times
-	ivec times;
-	//!Default constructor
-	str ( ivec ids0, ivec times0 ) : ids ( ids0 ), times ( times0 ) {
-		it_assert_debug ( times0.length() == ids0.length(), "Incompatible input" );
-	};
+class str
+{
+	public:
+		//! vector id ids (non-unique!)
+		ivec ids;
+		//! vector of times
+		ivec times;
+		//!Default constructor
+		str (ivec ids0, ivec times0) : ids (ids0), times (times0) {
+			it_assert_debug (times0.length() == ids0.length(), "Incompatible input");
+		};
 };
 
@@ -82,148 +84,149 @@
 */
 
-class RV : public root {
-protected:
-	//! size of the data vector
-	int dsize;
-	//! number of individual rvs
-	int len;
-	//! Vector of unique IDs
-	ivec ids;
-	//! Vector of shifts from current time
-	ivec times;
-
-private:
-	//! auxiliary function used in constructor
-	void init ( const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times );
-	int init ( const string &name, int size );
-public:
-	//! \name Constructors
-	//!@{
-
-	//! Full constructor
-	RV ( const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times ) {
-		init ( in_names, in_sizes, in_times );
-	}
-
-	//! Constructor with times=0
-	RV ( const Array<std::string> &in_names, const ivec &in_sizes ) {
-		init ( in_names, in_sizes, zeros_i ( in_names.length() ) );
-	}
-
-	//! Constructor with sizes=1, times=0
-	RV ( const Array<std::string> &in_names ) {
-		init ( in_names, ones_i ( in_names.length() ), zeros_i ( in_names.length() ) );
-	}
-
-	//! Constructor of empty RV
-	RV() : dsize ( 0 ), len ( 0 ), ids ( 0 ), times ( 0 ) {}
-	//! Constructor of a single RV with given id
-	RV ( string name, int sz, int tm = 0 );
-	//!@}
-
-	//! \name Access functions
-	//!@{
-
-	//! State output, e.g. for debugging.
-	friend std::ostream &operator<< ( std::ostream &os, const RV &rv );
-
-	int _dsize() const {
-		return dsize;
-	}
-
-	//! Recount size of the corresponding data vector
-	int countsize() const;
-	ivec cumsizes() const;
-	int length() const {
-		return len;
-	}
-	int id ( int at ) const {
-		return ids ( at );
-	}
-	int size ( int at ) const {
-		return RV_SIZES ( ids ( at ) );
-	}
-	int time ( int at ) const {
-		return times ( at );
-	}
-	std::string name ( int at ) const {
-		return RV_NAMES ( ids ( at ) );
-	}
-	void set_time ( int at, int time0 ) {
-		times ( at ) = time0;
-	}
-	//!@}
-
-	//TODO why not inline and later??
-
-	//! \name Algebra on Random Variables
-	//!@{
-
-	//! Find indices of self in another rv, \return ivec of the same size as self.
-	ivec findself ( const RV &rv2 ) const;
-	//! Compare if \c rv2 is identical to this \c RV
-	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 );
-	//! Subtract  another variable from the current one
-	RV subt ( const RV &rv2 ) const;
-	//! Select only variables at indices ind
-	RV subselect ( const ivec &ind ) const;
-
-	//! Select only variables at indices ind
-	RV operator() ( const ivec &ind ) const {
-		return subselect ( ind );
-	}
-
-	//! Select from data vector starting at di1 to di2
-	RV operator() ( int di1, int di2 ) const;
-
-	//! Shift \c time by delta.
-	void t ( int delta );
-	//!@}
-
-	//!\name Relation to vectors
-	//!@{
-
-	//! generate \c str from rv, by expanding sizes TODO to_string..
-	str tostr() const;
-	//! when this rv is a part of bigger rv, this function returns indices 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 indices when copying data between 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;
-	//! Minimum time-offset
-	int mint() const {
-		return min ( times );
-	}
-	//!@}
-
-	// TODO aktualizovat dle soucasneho UI
-	/*! \brief UI for class RV (description of data vectors)
-
-	\code
-	rv = {
-	  type = "rv"; //identifier of the description
-	  // UNIQUE IDENTIFIER same names = same variable
-	  names = ["a", "b", "c", ...];   // which will be used e.g. in loggers
-
-	  //optional arguments
-	  sizes = [1, 2, 3, ...];         // (optional) default = ones()
-	  times = [-1, -2, 0, ...];       // time shifts with respect to current time (optional) default = zeros()
-	}
-	\endcode
-	*/
-	void from_setting ( const Setting &set );
-
-	// TODO dodelat void to_setting( Setting &set ) const;
-
-	//! Invalidate all named RVs. Use before initializing any RV instances, with care...
-	static void clear_all();
-};
-UIREGISTER ( RV );
+class RV : public root
+{
+	protected:
+		//! size of the data vector
+		int dsize;
+		//! number of individual rvs
+		int len;
+		//! Vector of unique IDs
+		ivec ids;
+		//! Vector of shifts from current time
+		ivec times;
+
+	private:
+		//! auxiliary function used in constructor
+		void init (const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times);
+		int init (const string &name, int size);
+	public:
+		//! \name Constructors
+		//!@{
+
+		//! Full constructor
+		RV (const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times) {
+			init (in_names, in_sizes, in_times);
+		}
+
+		//! Constructor with times=0
+		RV (const Array<std::string> &in_names, const ivec &in_sizes) {
+			init (in_names, in_sizes, zeros_i (in_names.length()));
+		}
+
+		//! Constructor with sizes=1, times=0
+		RV (const Array<std::string> &in_names) {
+			init (in_names, ones_i (in_names.length()), zeros_i (in_names.length()));
+		}
+
+		//! Constructor of empty RV
+		RV() : dsize (0), len (0), ids (0), times (0) {}
+		//! Constructor of a single RV with given id
+		RV (string name, int sz, int tm = 0);
+		//!@}
+
+		//! \name Access functions
+		//!@{
+
+		//! State output, e.g. for debugging.
+		friend std::ostream &operator<< (std::ostream &os, const RV &rv);
+
+		int _dsize() const {
+			return dsize;
+		}
+
+		//! Recount size of the corresponding data vector
+		int countsize() const;
+		ivec cumsizes() const;
+		int length() const {
+			return len;
+		}
+		int id (int at) const {
+			return ids (at);
+		}
+		int size (int at) const {
+			return RV_SIZES (ids (at));
+		}
+		int time (int at) const {
+			return times (at);
+		}
+		std::string name (int at) const {
+			return RV_NAMES (ids (at));
+		}
+		void set_time (int at, int time0) {
+			times (at) = time0;
+		}
+		//!@}
+
+		//TODO why not inline and later??
+
+		//! \name Algebra on Random Variables
+		//!@{
+
+		//! Find indices of self in another rv, \return ivec of the same size as self.
+		ivec findself (const RV &rv2) const;
+		//! Compare if \c rv2 is identical to this \c RV
+		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);
+		//! Subtract  another variable from the current one
+		RV subt (const RV &rv2) const;
+		//! Select only variables at indices ind
+		RV subselect (const ivec &ind) const;
+
+		//! Select only variables at indices ind
+		RV operator() (const ivec &ind) const {
+			return subselect (ind);
+		}
+
+		//! Select from data vector starting at di1 to di2
+		RV operator() (int di1, int di2) const;
+
+		//! Shift \c time by delta.
+		void t (int delta);
+		//!@}
+
+		//!\name Relation to vectors
+		//!@{
+
+		//! generate \c str from rv, by expanding sizes TODO to_string..
+		str tostr() const;
+		//! when this rv is a part of bigger rv, this function returns indices 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 indices when copying data between 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;
+		//! Minimum time-offset
+		int mint() const {
+			return min (times);
+		}
+		//!@}
+
+		// TODO aktualizovat dle soucasneho UI
+		/*! \brief UI for class RV (description of data vectors)
+
+		\code
+		rv = {
+		  type = "rv"; //identifier of the description
+		  // UNIQUE IDENTIFIER same names = same variable
+		  names = ["a", "b", "c", ...];   // which will be used e.g. in loggers
+
+		  //optional arguments
+		  sizes = [1, 2, 3, ...];         // (optional) default = ones()
+		  times = [-1, -2, 0, ...];       // time shifts with respect to current time (optional) default = zeros()
+		}
+		\endcode
+		*/
+		void from_setting (const Setting &set);
+
+		// TODO dodelat void to_setting( Setting &set ) const;
+
+		//! Invalidate all named RVs. Use before initializing any RV instances, with care...
+		static void clear_all();
+};
+UIREGISTER (RV);
 
 //! Concat two random variables
-RV concat ( const RV &rv1, const RV &rv2 );
+RV concat (const RV &rv1, const RV &rv2);
 
 //!Default empty RV that can be used as default argument
@@ -232,23 +235,24 @@
 //! Class representing function \f$f(x)\f$ of variable \f$x\f$ represented by \c rv
 
-class fnc : public root {
-protected:
-	//! Length of the output vector
-	int dimy;
-public:
-	//!default constructor
-	fnc() {};
-	//! function evaluates numerical value of \f$f(x)\f$ at \f$x=\f$ \c cond
-	virtual vec eval ( const vec &cond ) {
-		return vec ( 0 );
-	};
-
-	//! function substitutes given value into an appropriate position
-	virtual void condition ( const vec &val ) {};
-
-	//! access function
-	int dimension() const {
-		return dimy;
-	}
+class fnc : public root
+{
+	protected:
+		//! Length of the output vector
+		int dimy;
+	public:
+		//!default constructor
+		fnc() {};
+		//! function evaluates numerical value of \f$f(x)\f$ at \f$x=\f$ \c cond
+		virtual vec eval (const vec &cond) {
+			return vec (0);
+		};
+
+		//! function substitutes given value into an appropriate position
+		virtual void condition (const vec &val) {};
+
+		//! access function
+		int dimension() const {
+			return dimy;
+		}
 };
 
@@ -257,140 +261,141 @@
 //! Probability density function with numerical statistics, e.g. posterior density.
 
-class epdf : public root {
-protected:
-	//! dimension of the random variable
-	int dim;
-	//! Description of the random variable
-	RV rv;
-
-public:
-	/*! \name Constructors
-	 Construction of each epdf should support two types of constructors:
-	\li empty constructor,
-	\li copy constructor,
-
-	The following constructors should be supported for convenience:
-	\li constructor followed by calling \c set_parameters()
-	\li constructor accepting random variables calling \c set_rv()
-
-	 All internal data structures are constructed as empty. Their values (including sizes) will be set by method \c set_parameters(). This way references can be initialized in constructors.
-	@{*/
-	epdf() : dim ( 0 ), rv() {};
-	epdf ( const epdf &e ) : dim ( e.dim ), rv ( e.rv ) {};
-	epdf ( const RV &rv0 ) : dim ( rv0._dsize() ) {
-		set_rv ( rv0 );
-	};
-	void set_parameters ( int dim0 ) {
-		dim = dim0;
-	}
-	//!@}
-
-	//! \name Matematical Operations
-	//!@{
-
-	//! Returns a sample, \f$ x \f$ from density \f$ f_x()\f$
-	virtual vec sample() const {
-		it_error ( "not implemneted" );
-		return vec ( 0 );
-	};
-	//! Returns N samples, \f$ [x_1 , x_2 , \ldots \ \f$  from density \f$ f_x(rv)\f$
-	virtual mat sample_m ( int N ) const;
-	//! Compute log-probability of argument \c val
-	//! In case the argument is out of suport return -Infinity
-	virtual double evallog ( const vec &val ) const {
-		it_error ( "not implemneted" );
-		return 0.0;
-	};
-	//! 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;
-	}
-	//! 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
-	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;
-	}
-
-	//! return expected value
-	virtual vec mean() const {
-		it_error ( "not implemneted" );
-		return vec ( 0 );
-	};
-
-	//! return expected variance (not covariance!)
-	virtual vec variance() const {
-		it_error ( "not implemneted" );
-		return vec ( 0 );
-	};
-	//! Lower and upper bounds of \c percentage % quantile, returns mean-2*sigma as default
-	virtual void qbounds ( vec &lb, vec &ub, double percentage = 0.95 ) const {
-		vec mea = mean();
-		vec std = sqrt ( variance() );
-		lb = mea - 2 * std;
-		ub = mea + 2 * std;
-	};
-	//!@}
-
-	//! \name Connection to other classes
-	//! Description of the random quantity via attribute \c rv is optional.
-	//! For operations such as sampling \c rv does not need to be set. However, for \c marginalization
-	//! and \c conditioning \c rv has to be set. NB:
-	//! @{
-
-	//!Name its rv
-	void set_rv ( const RV &rv0 ) {
-		rv = rv0;
-	}   //it_assert_debug(isnamed(),""); };
-	//! True if rv is assigned
-	bool isnamed() const {
-		bool b = ( dim == rv._dsize() );
-		return b;
-	}
-	//! Return name (fails when isnamed is false)
-	const RV& _rv() const {
-		it_assert_debug ( isnamed(), "" );
-		return rv;
-	}
-	//!@}
-
-	//! \name Access to attributes
-	//! @{
-
-	//! 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 ) {
-		RV* r = UI::build<RV> ( set, "rv" );
-		if ( r ) {
-			set_rv ( *r );
-			delete r;
-		}
-	}
+class epdf : public root
+{
+	protected:
+		//! dimension of the random variable
+		int dim;
+		//! Description of the random variable
+		RV rv;
+
+	public:
+		/*! \name Constructors
+		 Construction of each epdf should support two types of constructors:
+		\li empty constructor,
+		\li copy constructor,
+
+		The following constructors should be supported for convenience:
+		\li constructor followed by calling \c set_parameters()
+		\li constructor accepting random variables calling \c set_rv()
+
+		 All internal data structures are constructed as empty. Their values (including sizes) will be set by method \c set_parameters(). This way references can be initialized in constructors.
+		@{*/
+		epdf() : dim (0), rv() {};
+		epdf (const epdf &e) : dim (e.dim), rv (e.rv) {};
+		epdf (const RV &rv0) : dim (rv0._dsize()) {
+			set_rv (rv0);
+		};
+		void set_parameters (int dim0) {
+			dim = dim0;
+		}
+		//!@}
+
+		//! \name Matematical Operations
+		//!@{
+
+		//! Returns a sample, \f$ x \f$ from density \f$ f_x()\f$
+		virtual vec sample() const {
+			it_error ("not implemneted");
+			return vec (0);
+		};
+		//! Returns N samples, \f$ [x_1 , x_2 , \ldots \ \f$  from density \f$ f_x(rv)\f$
+		virtual mat sample_m (int N) const;
+		//! Compute log-probability of argument \c val
+		//! In case the argument is out of suport return -Infinity
+		virtual double evallog (const vec &val) const {
+			it_error ("not implemneted");
+			return 0.0;
+		};
+		//! 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;
+		}
+		//! 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
+		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;
+		}
+
+		//! return expected value
+		virtual vec mean() const {
+			it_error ("not implemneted");
+			return vec (0);
+		};
+
+		//! return expected variance (not covariance!)
+		virtual vec variance() const {
+			it_error ("not implemneted");
+			return vec (0);
+		};
+		//! Lower and upper bounds of \c percentage % quantile, returns mean-2*sigma as default
+		virtual void qbounds (vec &lb, vec &ub, double percentage = 0.95) const {
+			vec mea = mean();
+			vec std = sqrt (variance());
+			lb = mea - 2 * std;
+			ub = mea + 2 * std;
+		};
+		//!@}
+
+		//! \name Connection to other classes
+		//! Description of the random quantity via attribute \c rv is optional.
+		//! For operations such as sampling \c rv does not need to be set. However, for \c marginalization
+		//! and \c conditioning \c rv has to be set. NB:
+		//! @{
+
+		//!Name its rv
+		void set_rv (const RV &rv0) {
+			rv = rv0;
+		}   //it_assert_debug(isnamed(),""); };
+		//! True if rv is assigned
+		bool isnamed() const {
+			bool b = (dim == rv._dsize());
+			return b;
+		}
+		//! Return name (fails when isnamed is false)
+		const RV& _rv() const {
+			it_assert_debug (isnamed(), "");
+			return rv;
+		}
+		//!@}
+
+		//! \name Access to attributes
+		//! @{
+
+		//! 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) {
+			RV* r = UI::build<RV> (set, "rv");
+			if (r) {
+				set_rv (*r);
+				delete r;
+			}
+		}
 
 };
@@ -399,115 +404,117 @@
 //! Conditional probability density, e.g. modeling some dependencies.
 //TODO Samplecond can be generalized
-class mpdf : public root {
-protected:
-	//!dimension of the condition
-	int dimc;
-	//! random variable in condition
-	RV rvc;
-private:
-	//! internal epdf, used only as cache to avoid virtual calls of \c _rv() and \c _dimension()
-	epdf* ep;
-	
-protected:
-	void set_ep(epdf &iepdf) {
-		ep = &iepdf;
-	}
-
-public:
-	//! \name Constructors
-	//! @{
-
-	mpdf() : dimc ( 0 ), rvc(), ep(NULL) { }
-
-	mpdf ( const mpdf &m ) : dimc ( m.dimc ), rvc ( m.rvc ), ep ( m.ep ) { }
-	//!@}
-
-	//! \name Matematical operations
-	//!@{
-
-	//! 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
-	virtual vec samplecond ( const vec &cond ){it_error("Not implemented");return vec(0);};
-
-	//! 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
-	virtual mat samplecond_m ( const vec &cond, int N ){it_error("Not implemented");return mat();}
-
-
-	//! Shortcut for conditioning and evaluation of the internal epdf. In some cases,  this operation can be implemented efficiently.
-	virtual double evallogcond ( const vec &dt, const vec &cond ){it_error("Not implemented");return 0.0;}
-
-	//! Matrix version of evallogcond
-	virtual vec evallogcond_m ( const mat &Dt, const vec &cond ){it_error("Not implemented");return vec(0);}
-
-	//! Array<vec> version of evallogcond
-	virtual vec evallogcond_m ( const Array<vec> &Dt, const vec &cond ){it_error("Not implemented");return vec(0);}
-
-	//! \name Access to attributes
-	//! @{
-
-	RV _rv() {
-		return ep->_rv();
-	}
-	RV _rvc() {
-		it_assert_debug ( isnamed(), "" );
-		return rvc;
-	}
-	int dimension() {
-		return ep->dimension();
-	}
-	int dimensionc() {
-		return dimc;
-	}
-
-	//! Load from structure with elements:
-	//!  \code
-	//! { class = "mpdf_offspring",
-	//!   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 );
-	//!@}
-
-	//! \name Connection to other objects
-	//!@{
-	void set_rvc ( const RV &rvc0 ) {
-		rvc = rvc0;
-	}
-	void set_rv ( const RV &rv0 ) {
-		ep->set_rv ( rv0 );
-	}
-	bool isnamed() {
-		return ( ep->isnamed() ) && ( dimc == rvc._dsize() );
-	}
-	//!@}
+class mpdf : public root
+{
+	protected:
+		//!dimension of the condition
+		int dimc;
+		//! random variable in condition
+		RV rvc;
+	private:
+		//! internal epdf, used only as cache to avoid virtual calls of \c _rv() and \c _dimension()
+		epdf* ep;
+
+	protected:
+		void set_ep (epdf &iepdf) {
+			ep = &iepdf;
+		}
+
+	public:
+		//! \name Constructors
+		//! @{
+
+		mpdf() : dimc (0), rvc(), ep (NULL) { }
+
+		mpdf (const mpdf &m) : dimc (m.dimc), rvc (m.rvc), ep (m.ep) { }
+		//!@}
+
+		//! \name Matematical operations
+		//!@{
+
+		//! 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
+		virtual vec samplecond (const vec &cond) {it_error ("Not implemented");return vec (0);};
+
+		//! 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
+		virtual mat samplecond_m (const vec &cond, int N) {it_error ("Not implemented");return mat();}
+
+
+		//! Shortcut for conditioning and evaluation of the internal epdf. In some cases,  this operation can be implemented efficiently.
+		virtual double evallogcond (const vec &dt, const vec &cond) {it_error ("Not implemented");return 0.0;}
+
+		//! Matrix version of evallogcond
+		virtual vec evallogcond_m (const mat &Dt, const vec &cond) {it_error ("Not implemented");return vec (0);}
+
+		//! Array<vec> version of evallogcond
+		virtual vec evallogcond_m (const Array<vec> &Dt, const vec &cond) {it_error ("Not implemented");return vec (0);}
+
+		//! \name Access to attributes
+		//! @{
+
+		RV _rv() {
+			return ep->_rv();
+		}
+		RV _rvc() {
+			it_assert_debug (isnamed(), "");
+			return rvc;
+		}
+		int dimension() {
+			return ep->dimension();
+		}
+		int dimensionc() {
+			return dimc;
+		}
+
+		//! Load from structure with elements:
+		//!  \code
+		//! { class = "mpdf_offspring",
+		//!   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);
+		//!@}
+
+		//! \name Connection to other objects
+		//!@{
+		void set_rvc (const RV &rvc0) {
+			rvc = rvc0;
+		}
+		void set_rv (const RV &rv0) {
+			ep->set_rv (rv0);
+		}
+		bool isnamed() {
+			return (ep->isnamed()) && (dimc == rvc._dsize());
+		}
+		//!@}
 };
 
 template <class EPDF>
-class mpdf_internal: public mpdf{
+class mpdf_internal: public mpdf
+{
 	protected :
 		EPDF iepdf;
 	public:
 		//! constructor
-		mpdf_internal(): mpdf(),iepdf(){set_ep(iepdf);}
+		mpdf_internal() : mpdf(), iepdf() {set_ep (iepdf);}
 		//! Update \c iepdf so that it represents this mpdf conditioned on \c rvc = cond
 		//! This function provides convenient reimplementation in offsprings
-			virtual void condition ( const vec &cond ) {
-				it_error ( "Not implemented" );
-			};
+		virtual void condition (const vec &cond) {
+			it_error ("Not implemented");
+		};
 		//!access function to iepdf
-		EPDF& e(){return iepdf;}
-			
+		EPDF& e() {return iepdf;}
+
 		//! Reimplements samplecond using \c condition()
-		vec samplecond ( const vec &cond );
+		vec samplecond (const vec &cond);
 		//! Reimplements evallogcond using \c condition()
-		double evallogcond ( const vec &val, const vec &cond );
-		//! Efficient version of evallogcond for matrices 
-		virtual vec evallogcond_m ( const mat &Dt, const vec &cond );
+		double evallogcond (const vec &val, const vec &cond);
+		//! Efficient version of evallogcond for matrices
+		virtual vec evallogcond_m (const mat &Dt, const vec &cond);
 		//! Efficient version of evallogcond for Array<vec>
-		virtual vec evallogcond_m ( const Array<vec> &Dt, const vec &cond );
-		//! Efficient version of samplecond 
-		virtual mat samplecond_m ( const vec &cond, int N );
+		virtual vec evallogcond_m (const Array<vec> &Dt, const vec &cond);
+		//! Efficient version of samplecond
+		virtual mat samplecond_m (const vec &cond, int N);
 };
 
@@ -537,119 +544,122 @@
 
 */
-class datalink {
-protected:
-	//! Remember how long val should be
-	int downsize;
-
-	//! Remember how long val of "Up" should be
-	int upsize;
-
-	//! val-to-val link, indices of the upper val
-	ivec v2v_up;
-
-public:
-	//! Constructor
-	datalink() : downsize ( 0 ), upsize ( 0 ) { }
-	datalink ( const RV &rv, const RV &rv_up ) {
-		set_connection ( rv, rv_up );
-	}
-
-	//! set connection, rv must be fully present in rv_up
-	void set_connection ( const RV &rv, const RV &rv_up ) {
-		downsize = rv._dsize();
-		upsize = rv_up._dsize();
-		v2v_up = rv.dataind ( rv_up );
-
-		it_assert_debug ( v2v_up.length() == downsize, "rv is not fully in rv_up" );
-	}
-
-	//! set connection using indices
-	void set_connection ( int ds, int us, const ivec &upind ) {
-		downsize = ds;
-		upsize = us;
-		v2v_up = upind;
-
-		it_assert_debug ( v2v_up.length() == downsize, "rv is not fully in rv_up" );
-	}
-
-	//! Get val for myself from val of "Up"
-	vec pushdown ( const vec &val_up ) {
-		it_assert_debug ( upsize == val_up.length(), "Wrong val_up" );
-		return get_vec ( val_up, v2v_up );
-	}
-
-	//! Fill val of "Up" by my pieces
-	void pushup ( vec &val_up, const vec &val ) {
-		it_assert_debug ( downsize == val.length(), "Wrong val" );
-		it_assert_debug ( upsize == val_up.length(), "Wrong val_up" );
-		set_subvector ( val_up, v2v_up, val );
-	}
+class datalink
+{
+	protected:
+		//! Remember how long val should be
+		int downsize;
+
+		//! Remember how long val of "Up" should be
+		int upsize;
+
+		//! val-to-val link, indices of the upper val
+		ivec v2v_up;
+
+	public:
+		//! Constructor
+		datalink() : downsize (0), upsize (0) { }
+		datalink (const RV &rv, const RV &rv_up) {
+			set_connection (rv, rv_up);
+		}
+
+		//! set connection, rv must be fully present in rv_up
+		void set_connection (const RV &rv, const RV &rv_up) {
+			downsize = rv._dsize();
+			upsize = rv_up._dsize();
+			v2v_up = rv.dataind (rv_up);
+
+			it_assert_debug (v2v_up.length() == downsize, "rv is not fully in rv_up");
+		}
+
+		//! set connection using indices
+		void set_connection (int ds, int us, const ivec &upind) {
+			downsize = ds;
+			upsize = us;
+			v2v_up = upind;
+
+			it_assert_debug (v2v_up.length() == downsize, "rv is not fully in rv_up");
+		}
+
+		//! Get val for myself from val of "Up"
+		vec pushdown (const vec &val_up) {
+			it_assert_debug (upsize == val_up.length(), "Wrong val_up");
+			return get_vec (val_up, v2v_up);
+		}
+
+		//! Fill val of "Up" by my pieces
+		void pushup (vec &val_up, const vec &val) {
+			it_assert_debug (downsize == val.length(), "Wrong val");
+			it_assert_debug (upsize == val_up.length(), "Wrong val_up");
+			set_subvector (val_up, v2v_up, val);
+		}
 };
 
 //! Data link with a condition.
-class datalink_m2e: public datalink {
-protected:
-	//! Remember how long cond should be
-	int condsize;
-
-	//!upper_val-to-local_cond link, indices of the upper val
-	ivec v2c_up;
-
-	//!upper_val-to-local_cond link, indices of the local cond
-	ivec v2c_lo;
-
-public:
-	//! Constructor
-	datalink_m2e() : condsize ( 0 ) { }
-
-	void set_connection ( const RV &rv, const RV &rvc, const RV &rv_up ) {
-		datalink::set_connection ( rv, rv_up );
-		condsize = rvc._dsize();
-		//establish v2c connection
-		rvc.dataind ( rv_up, v2c_lo, v2c_up );
-	}
-
-	//!Construct condition
-	vec get_cond ( const vec &val_up ) {
-		vec tmp ( condsize );
-		set_subvector ( tmp, v2c_lo, val_up ( v2c_up ) );
-		return tmp;
-	}
-
-	void pushup_cond ( vec &val_up, const vec &val, const vec &cond ) {
-		it_assert_debug ( downsize == val.length(), "Wrong val" );
-		it_assert_debug ( upsize == val_up.length(), "Wrong val_up" );
-		set_subvector ( val_up, v2v_up, val );
-		set_subvector ( val_up, v2c_up, cond );
-	}
+class datalink_m2e: public datalink
+{
+	protected:
+		//! Remember how long cond should be
+		int condsize;
+
+		//!upper_val-to-local_cond link, indices of the upper val
+		ivec v2c_up;
+
+		//!upper_val-to-local_cond link, indices of the local cond
+		ivec v2c_lo;
+
+	public:
+		//! Constructor
+		datalink_m2e() : condsize (0) { }
+
+		void set_connection (const RV &rv, const RV &rvc, const RV &rv_up) {
+			datalink::set_connection (rv, rv_up);
+			condsize = rvc._dsize();
+			//establish v2c connection
+			rvc.dataind (rv_up, v2c_lo, v2c_up);
+		}
+
+		//!Construct condition
+		vec get_cond (const vec &val_up) {
+			vec tmp (condsize);
+			set_subvector (tmp, v2c_lo, val_up (v2c_up));
+			return tmp;
+		}
+
+		void pushup_cond (vec &val_up, const vec &val, const vec &cond) {
+			it_assert_debug (downsize == val.length(), "Wrong val");
+			it_assert_debug (upsize == val_up.length(), "Wrong val_up");
+			set_subvector (val_up, v2v_up, val);
+			set_subvector (val_up, v2c_up, cond);
+		}
 };
 
 //!DataLink is a connection between mpdf and its superordinate (Up)
 //! This class links
-class datalink_m2m: public datalink_m2e {
-protected:
-	//!cond-to-cond link, indices of the upper cond
-	ivec c2c_up;
-	//!cond-to-cond link, indices of the local cond
-	ivec c2c_lo;
-
-public:
-	//! Constructor
-	datalink_m2m() {};
-	void set_connection ( const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up ) {
-		datalink_m2e::set_connection ( rv, rvc, rv_up );
-		//establish c2c connection
-		rvc.dataind ( rvc_up, c2c_lo, c2c_up );
-		it_assert_debug ( c2c_lo.length() + v2c_lo.length() == condsize, "cond is not fully given" );
-	}
-
-	//! Get cond for myself from val and cond of "Up"
-	vec get_cond ( const vec &val_up, const vec &cond_up ) {
-		vec tmp ( condsize );
-		set_subvector ( tmp, v2c_lo, val_up ( v2c_up ) );
-		set_subvector ( tmp, c2c_lo, cond_up ( c2c_up ) );
-		return tmp;
-	}
-	//! Fill
+class datalink_m2m: public datalink_m2e
+{
+	protected:
+		//!cond-to-cond link, indices of the upper cond
+		ivec c2c_up;
+		//!cond-to-cond link, indices of the local cond
+		ivec c2c_lo;
+
+	public:
+		//! Constructor
+		datalink_m2m() {};
+		void set_connection (const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up) {
+			datalink_m2e::set_connection (rv, rvc, rv_up);
+			//establish c2c connection
+			rvc.dataind (rvc_up, c2c_lo, c2c_up);
+			it_assert_debug (c2c_lo.length() + v2c_lo.length() == condsize, "cond is not fully given");
+		}
+
+		//! Get cond for myself from val and cond of "Up"
+		vec get_cond (const vec &val_up, const vec &cond_up) {
+			vec tmp (condsize);
+			set_subvector (tmp, v2c_lo, val_up (v2c_up));
+			set_subvector (tmp, c2c_lo, cond_up (c2c_up));
+			return tmp;
+		}
+		//! Fill
 
 };
@@ -660,42 +670,43 @@
 This class abstracts logging of results from implementation. This class replaces direct logging of results (e.g. to files or to global variables) by calling methods of a logger. Specializations of this abstract class for specific storage method are designed.
  */
-class logger : public root {
-protected:
-	//! RVs of all logged variables.
-	Array<RV> entries;
-	//! Names of logged quantities, e.g. names of algorithm variants
-	Array<string> names;
-public:
-	//!Default constructor
-	logger() : entries ( 0 ), names ( 0 ) {}
-
-	//! returns an identifier which will be later needed for calling the \c logit() function
-	//! For empty RV it returns -1, this entry will be ignored by \c logit().
-	virtual int add ( const RV &rv, string prefix = "" ) {
-		int id;
-		if ( rv._dsize() > 0 ) {
-			id = entries.length();
-			names = concat ( names, prefix ); // diff
-			entries.set_length ( id + 1, true );
-			entries ( id ) = rv;
-		} else {
-			id = -1;
-		}
-		return id; // identifier of the last entry
-	}
-
-	//! log this vector
-	virtual void logit ( int id, const vec &v ) = 0;
-	//! log this double
-	virtual void logit ( int id, const double &d ) = 0;
-
-	//! Shifts storage position for another time step.
-	virtual void step() = 0;
-
-	//! Finalize storing information
-	virtual void finalize() {};
-
-	//! Initialize the storage
-	virtual void init() {};
+class logger : public root
+{
+	protected:
+		//! RVs of all logged variables.
+		Array<RV> entries;
+		//! Names of logged quantities, e.g. names of algorithm variants
+		Array<string> names;
+	public:
+		//!Default constructor
+		logger() : entries (0), names (0) {}
+
+		//! returns an identifier which will be later needed for calling the \c logit() function
+		//! For empty RV it returns -1, this entry will be ignored by \c logit().
+		virtual int add (const RV &rv, string prefix = "") {
+			int id;
+			if (rv._dsize() > 0) {
+				id = entries.length();
+				names = concat (names, prefix);   // diff
+				entries.set_length (id + 1, true);
+				entries (id) = rv;
+			} else {
+				id = -1;
+			}
+			return id; // identifier of the last entry
+		}
+
+		//! log this vector
+		virtual void logit (int id, const vec &v) = 0;
+		//! log this double
+		virtual void logit (int id, const double &d) = 0;
+
+		//! Shifts storage position for another time step.
+		virtual void step() = 0;
+
+		//! Finalize storing information
+		virtual void finalize() {};
+
+		//! Initialize the storage
+		virtual void init() {};
 
 };
@@ -704,57 +715,59 @@
 
 */
-class mepdf : public mpdf {
-
-shared_ptr<epdf> ipdf;
-public:
-	//!Default constructor
-	mepdf() { }
-
-	mepdf ( shared_ptr<epdf> em ) {
-		ipdf = em;
-		set_ep(*ipdf.get());
-		dimc = 0;
-	}
-
-	//! empty
-	void condition ( const vec &cond );
-
-	//! Load from structure with elements:
-	//!  \code
-	//! { class = "mepdf",
-	//!   epdf = {class="epdf_offspring",...}
-	//! }
-	//! \endcode
-	//!@}
-	void from_setting ( const Setting &set );
-};
-UIREGISTER ( mepdf );
+class mepdf : public mpdf
+{
+
+		shared_ptr<epdf> ipdf;
+	public:
+		//!Default constructor
+		mepdf() { }
+
+		mepdf (shared_ptr<epdf> em) {
+			ipdf = em;
+			set_ep (*ipdf.get());
+			dimc = 0;
+		}
+
+		//! empty
+		void condition (const vec &cond);
+
+		//! Load from structure with elements:
+		//!  \code
+		//! { class = "mepdf",
+		//!   epdf = {class="epdf_offspring",...}
+		//! }
+		//! \endcode
+		//!@}
+		void from_setting (const Setting &set);
+};
+UIREGISTER (mepdf);
 
 //!\brief Chain rule of pdfs - abstract part common for mprod and merger.
 //!this abstract class is common to epdf and mpdf
 //!\todo Think of better design - global functions rv=get_rv(Array<mpdf*> mpdfs); ??
-class compositepdf {
-protected:
-	//! Elements of composition
-	Array<mpdf*> mpdfs;
-	bool owning_mpdfs;
-public:
-	compositepdf() : mpdfs ( 0 ) {};
-	compositepdf ( Array<mpdf*> A0, bool own = false ) {
-		set_elements ( A0, own );
-	};
-	void set_elements ( Array<mpdf*> A0, bool own = false ) {
-		mpdfs = A0;
-		owning_mpdfs = own;
-	};
-	//! 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 );
-	~compositepdf() {
-		if ( owning_mpdfs ) for ( int i = 0; i < mpdfs.length(); i++ ) {
-				delete mpdfs ( i );
-			}
-	};
+class compositepdf
+{
+	protected:
+		//! Elements of composition
+		Array<mpdf*> mpdfs;
+		bool owning_mpdfs;
+	public:
+		compositepdf() : mpdfs (0) {};
+		compositepdf (Array<mpdf*> A0, bool own = false) {
+			set_elements (A0, own);
+		};
+		void set_elements (Array<mpdf*> A0, bool own = false) {
+			mpdfs = A0;
+			owning_mpdfs = own;
+		};
+		//! 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);
+		~compositepdf() {
+			if (owning_mpdfs) for (int i = 0; i < mpdfs.length(); i++) {
+					delete mpdfs (i);
+				}
+		};
 };
 
@@ -766,67 +779,68 @@
 */
 
-class DS : public root {
-protected:
-	int dtsize;
-	int utsize;
-	//!Description of data returned by \c getdata().
-	RV Drv;
-	//!Description of data witten by by \c write().
-	RV Urv; //
-	//! Remember its own index in Logger L
-	int L_dt, L_ut;
-public:
-	//! default constructors
-	DS() : Drv(), Urv() {};
-	//! Returns full vector of observed data=[output, input]
-	virtual void getdata ( vec &dt ) {
-		it_error ( "abstract class" );
-	};
-	//! Returns data records at indeces.
-	virtual void getdata ( vec &dt, const ivec &indeces ) {
-		it_error ( "abstract class" );
-	};
-	//! Accepts action variable and schedule it for application.
-	virtual void write ( vec &ut ) {
-		it_error ( "abstract class" );
-	};
-	//! Accepts action variables at specific indeces
-	virtual void write ( vec &ut, const ivec &indeces ) {
-		it_error ( "abstract class" );
-	};
-
-	//! Moves from \f$ t \f$ to \f$ t+1 \f$, i.e. perfroms the actions and reads response of the system.
-	virtual void step() = 0;
-
-	//! Register DS for logging into logger L
-	virtual void log_add ( logger &L ) {
-		it_assert_debug ( dtsize == Drv._dsize(), "" );
-		it_assert_debug ( utsize == Urv._dsize(), "" );
-
-		L_dt = L.add ( Drv, "" );
-		L_ut = L.add ( Urv, "" );
-	}
-	//! Register DS for logging into logger L
-	virtual void logit ( logger &L ) {
-		vec tmp ( Drv._dsize() + Urv._dsize() );
-		getdata ( tmp );
-		// d is first in getdata
-		L.logit ( L_dt, tmp.left ( Drv._dsize() ) );
-		// u follows after d in getdata
-		L.logit ( L_ut, tmp.mid ( Drv._dsize(), Urv._dsize() ) );
-	}
-	//!access function
-	virtual RV _drv() const {
-		return concat ( Drv, Urv );
-	}
-	//!access function
-	const RV& _urv() const {
-		return Urv;
-	}
-	//! set random rvariables
-	virtual void set_drv ( const  RV &drv, const RV &urv ) {
-		Drv = drv;
-		Urv = urv;
-	}
+class DS : public root
+{
+	protected:
+		int dtsize;
+		int utsize;
+		//!Description of data returned by \c getdata().
+		RV Drv;
+		//!Description of data witten by by \c write().
+		RV Urv; //
+		//! Remember its own index in Logger L
+		int L_dt, L_ut;
+	public:
+		//! default constructors
+		DS() : Drv(), Urv() {};
+		//! Returns full vector of observed data=[output, input]
+		virtual void getdata (vec &dt) {
+			it_error ("abstract class");
+		};
+		//! Returns data records at indeces.
+		virtual void getdata (vec &dt, const ivec &indeces) {
+			it_error ("abstract class");
+		};
+		//! Accepts action variable and schedule it for application.
+		virtual void write (vec &ut) {
+			it_error ("abstract class");
+		};
+		//! Accepts action variables at specific indeces
+		virtual void write (vec &ut, const ivec &indeces) {
+			it_error ("abstract class");
+		};
+
+		//! Moves from \f$ t \f$ to \f$ t+1 \f$, i.e. perfroms the actions and reads response of the system.
+		virtual void step() = 0;
+
+		//! Register DS for logging into logger L
+		virtual void log_add (logger &L) {
+			it_assert_debug (dtsize == Drv._dsize(), "");
+			it_assert_debug (utsize == Urv._dsize(), "");
+
+			L_dt = L.add (Drv, "");
+			L_ut = L.add (Urv, "");
+		}
+		//! Register DS for logging into logger L
+		virtual void logit (logger &L) {
+			vec tmp (Drv._dsize() + Urv._dsize());
+			getdata (tmp);
+			// d is first in getdata
+			L.logit (L_dt, tmp.left (Drv._dsize()));
+			// u follows after d in getdata
+			L.logit (L_ut, tmp.mid (Drv._dsize(), Urv._dsize()));
+		}
+		//!access function
+		virtual RV _drv() const {
+			return concat (Drv, Urv);
+		}
+		//!access function
+		const RV& _urv() const {
+			return Urv;
+		}
+		//! set random rvariables
+		virtual void set_drv (const  RV &drv, const RV &urv) {
+			Drv = drv;
+			Urv = urv;
+		}
 };
 
@@ -852,155 +866,201 @@
 */
 
-class BM : public root {
-protected:
-	//! Random variable of the data (optional)
-	RV drv;
-	//!Logarithm of marginalized data likelihood.
-	double ll;
-	//!  If true, the filter will compute likelihood of the data record and store it in \c ll . Set to false if you want to save computational time.
-	bool evalll;
-public:
-	//! \name Constructors
-	//! @{
-
-	BM() : ll ( 0 ), evalll ( true ), LIDs ( 4 ), LFlags ( 4 ) {
-		LIDs = -1;/*empty IDs*/
-		LFlags = 0;
-		LFlags ( 0 ) = 1;/*log only mean*/
-	};
-	BM ( const BM &B ) :  drv ( B.drv ), ll ( B.ll ), evalll ( B.evalll ) {}
-	//! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually!
-	//! Prototype: \code BM* _copy_() const {return new BM(*this);} \endcode
-	virtual BM* _copy_() const {
-		return NULL;
-	};
-	//!@}
-
-	//! \name Mathematical operations
-	//!@{
-
-	/*! \brief Incremental Bayes rule
-	@param dt vector of input data
-	*/
-	virtual void bayes ( const vec &dt ) = 0;
-	//! Batch Bayes rule (columns of Dt are observations)
-	virtual void bayesB ( const mat &Dt );
-	//! 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;
+class BM : public root
+{
+	protected:
+		//! Random variable of the data (optional)
+		RV drv;
+		//!Logarithm of marginalized data likelihood.
+		double ll;
+		//!  If true, the filter will compute likelihood of the data record and store it in \c ll . Set to false if you want to save computational time.
+		bool evalll;
+	public:
+		//! \name Constructors
+		//! @{
+
+		BM() : ll (0), evalll (true), LIDs (4), LFlags (4) {
+			LIDs = -1;/*empty IDs*/
+			LFlags = 0;
+			LFlags (0) = 1;  /*log only mean*/
+		};
+		BM (const BM &B) :  drv (B.drv), ll (B.ll), evalll (B.evalll) {}
+		//! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually!
+		//! Prototype: \code BM* _copy_() const {return new BM(*this);} \endcode
+		virtual BM* _copy_() const {
+			return NULL;
+		};
+		//!@}
+
+		//! \name Mathematical operations
+		//!@{
+
+		/*! \brief Incremental Bayes rule
+		@param dt vector of input data
+		*/
+		virtual void bayes (const vec &dt) = 0;
+		//! Batch Bayes rule (columns of Dt are observations)
+		virtual void bayesB (const mat &Dt);
+		//! 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;
+		}
+		//! 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;
+		}
+
+		//!Constructs a predictive density \f$ f(d_{t+1} |d_{t}, \ldots d_{0}) \f$
+		virtual epdf* epredictor() const {
+			it_error ("Not implemented");
+			return NULL;
+		};
+		//!Constructs a conditional density 1-step ahead predictor \f$ f(d_{t+1} |d_{t+h-1}, \ldots d_{t})
+		virtual mpdf* predictor() const {
+			it_error ("Not implemented");
+			return NULL;
+		};
+		//!@}
+
+		//! \name Extension to conditional BM
+		//! This extension is useful e.g. in Marginalized Particle Filter (\ref bdm::MPF).
+		//! Alternatively, it can be used for automated connection to DS when the condition is observed
+		//!@{
+
+		//! Name of extension variable
+		RV rvc;
+		//! access function
+		const RV& _rvc() const {
+			return rvc;
+		}
+
+		//! Substitute \c val for \c rvc.
+		virtual void condition (const vec &val) {
+			it_error ("Not implemented!");
+		};
+
+		//!@}
+
+
+		//! \name Access to attributes
+		//!@{
+
+		const RV& _drv() const {
+			return drv;
+		}
+		void set_drv (const RV &rv) {
+			drv = rv;
+		}
+		void set_rv (const RV &rv) {
+			const_cast<epdf&> (posterior()).set_rv (rv);
+		}
+		double _ll() const {
+			return ll;
+		}
+		void set_evalll (bool evl0) {
+			evalll = evl0;
+		}
+		virtual const epdf& posterior() const = 0;
+		virtual const epdf* _e() const = 0;
+		//!@}
+
+		//! \name Logging of results
+		//!@{
+
+		//! Set boolean options from a string, recognized are: "logbounds,logll"
+		virtual void set_options (const string &opt) {
+			LFlags (0) = 1;
+			if (opt.find ("logbounds") != string::npos) {
+				LFlags (1) = 1;
+				LFlags (2) = 1;
+			}
+			if (opt.find ("logll") != string::npos) {
+				LFlags (3) = 1;
+			}
+		}
+		//! IDs of storages in loggers 4:[1=mean,2=lb,3=ub,4=ll]
+		ivec LIDs;
+
+		//! Flags for logging - same size as LIDs, each entry correspond to the same in LIDs
+		ivec LFlags;
+		//! Add all logged variables to a logger
+		virtual void log_add (logger &L, const string &name = "") {
+			// internal
+			RV r;
+			if (posterior().isnamed()) {
+				r = posterior()._rv();
+			} else {
+				r = RV ("est", posterior().dimension());
+			};
+
+			// Add mean value
+			if (LFlags (0)) LIDs (0) = L.add (r, name + "mean_");
+			if (LFlags (1)) LIDs (1) = L.add (r, name + "lb_");
+			if (LFlags (2)) LIDs (2) = L.add (r, name + "ub_");
+			if (LFlags (3)) LIDs (3) = L.add (RV ("ll", 1), name);              //TODO: "local" RV
+		}
+		virtual void logit (logger &L) {
+			L.logit (LIDs (0), posterior().mean());
+			if (LFlags (1) || LFlags (2)) {        //if one of them is off, its LID==-1 and will not be stored
+				vec ub, lb;
+				posterior().qbounds (lb, ub);
+				L.logit (LIDs (1), lb);
+				L.logit (LIDs (2), ub);
+			}
+			if (LFlags (3)) L.logit (LIDs (3), ll);
+		}
+		//!@}
+};
+
+template<class EPDF>
+vec mpdf_internal<EPDF>::samplecond (const vec &cond)
+{
+	condition (cond);
+	vec temp = iepdf.sample();
+	return temp;
+}
+
+template<class EPDF>
+mat mpdf_internal<EPDF>::samplecond_m (const vec &cond, int N)
+{
+	condition (cond);
+	mat temp (dimension(), N);
+	vec smp (dimension());
+	for (int i = 0; i < N; i++) {
+		smp = iepdf.sample();
+		temp.set_col (i, smp);
 	}
-	//! 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;
-	}
-
-	//!Constructs a predictive density \f$ f(d_{t+1} |d_{t}, \ldots d_{0}) \f$
-	virtual epdf* epredictor() const {
-		it_error ( "Not implemented" );
-		return NULL;
-	};
-	//!Constructs a conditional density 1-step ahead predictor \f$ f(d_{t+1} |d_{t+h-1}, \ldots d_{t})
-	virtual mpdf* predictor() const {
-		it_error ( "Not implemented" );
-		return NULL;
-	};
-	//!@}
-
-	//! \name Extension to conditional BM
-	//! This extension is useful e.g. in Marginalized Particle Filter (\ref bdm::MPF).
-	//! Alternatively, it can be used for automated connection to DS when the condition is observed
-	//!@{
-
-	//! Name of extension variable
-	RV rvc;
-	//! access function
-	const RV& _rvc() const {
-		return rvc;
-	}
-
-	//! Substitute \c val for \c rvc.
-	virtual void condition ( const vec &val ) {
-		it_error ( "Not implemented!" );
-	};
-
-	//!@}
-
-
-	//! \name Access to attributes
-	//!@{
-
-	const RV& _drv() const {
-		return drv;
-	}
-	void set_drv ( const RV &rv ) {
-		drv = rv;
-	}
-	void set_rv ( const RV &rv ) {
-		const_cast<epdf&> ( posterior() ).set_rv ( rv );
-	}
-	double _ll() const {
-		return ll;
-	}
-	void set_evalll ( bool evl0 ) {
-		evalll = evl0;
-	}
-	virtual const epdf& posterior() const = 0;
-	virtual const epdf* _e() const = 0;
-	//!@}
-
-	//! \name Logging of results
-	//!@{
-
-	//! Set boolean options from a string, recognized are: "logbounds,logll"
-	virtual void set_options ( const string &opt ) {
-		LFlags ( 0 ) = 1;
-		if ( opt.find ( "logbounds" ) != string::npos ) {
-			LFlags ( 1 ) = 1;
-			LFlags ( 2 ) = 1;
-		}
-		if ( opt.find ( "logll" ) != string::npos ) {
-			LFlags ( 3 ) = 1;
-		}
-	}
-	//! IDs of storages in loggers 4:[1=mean,2=lb,3=ub,4=ll]
-	ivec LIDs;
-
-	//! Flags for logging - same size as LIDs, each entry correspond to the same in LIDs
-	ivec LFlags;
-	//! Add all logged variables to a logger
-	virtual void log_add ( logger &L, const string &name = "" ) {
-		// internal
-		RV r;
-		if ( posterior().isnamed() ) {
-			r = posterior()._rv();
-		} else {
-			r = RV ( "est", posterior().dimension() );
-		};
-
-		// Add mean value
-		if ( LFlags ( 0 ) ) LIDs ( 0 ) = L.add ( r, name + "mean_" );
-		if ( LFlags ( 1 ) ) LIDs ( 1 ) = L.add ( r, name + "lb_" );
-		if ( LFlags ( 2 ) ) LIDs ( 2 ) = L.add ( r, name + "ub_" );
-		if ( LFlags ( 3 ) ) LIDs ( 3 ) = L.add ( RV ( "ll", 1 ), name );    //TODO: "local" RV
-	}
-	virtual void logit ( logger &L ) {
-		L.logit ( LIDs ( 0 ), posterior().mean() );
-		if ( LFlags ( 1 ) || LFlags ( 2 ) ) {  //if one of them is off, its LID==-1 and will not be stored
-			vec ub, lb;
-			posterior().qbounds ( lb, ub );
-			L.logit ( LIDs ( 1 ), lb );
-			L.logit ( LIDs ( 2 ), ub );
-		}
-		if ( LFlags ( 3 ) ) L.logit ( LIDs ( 3 ), ll );
-	}
-	//!@}
-};
-
+
+	return temp;
+}
+
+template<class EPDF>
+double mpdf_internal<EPDF>::evallogcond (const vec &dt, const vec &cond)
+{
+	double tmp;
+	condition (cond);
+	tmp = iepdf.evallog (dt);
+	// it_assert_debug(std::isfinite(tmp), "Infinite value");
+	return tmp;
+}
+
+template<class EPDF>
+vec mpdf_internal<EPDF>::evallogcond_m (const mat &Dt, const vec &cond)
+{
+	condition (cond);
+	return iepdf.evallog_m (Dt);
+}
+
+template<class EPDF>
+vec mpdf_internal<EPDF>::evallogcond_m (const Array<vec> &Dt, const vec &cond)
+{
+	condition (cond);
+	return iepdf.evallog_m (Dt);
+}
 
 }; //namespace
Index: library/bdm/base/bdmbase.cpp
===================================================================
--- library/bdm/base/bdmbase.cpp (revision 487)
+++ library/bdm/base/bdmbase.cpp (revision 488)
@@ -140,44 +140,4 @@
 }
 
-template<class EPDF>
-vec mpdf_internal<EPDF>::samplecond ( const vec &cond ) {
-	condition ( cond );
-	vec temp = iepdf.sample();
-	return temp;
-}
-
-template<class EPDF>
-mat mpdf_internal<EPDF>::samplecond_m ( const vec &cond, int N ) {
-	condition ( cond );
-	mat temp ( dimension(), N );
-	vec smp ( dimension() );
-	for ( int i = 0; i < N; i++ ) {
-		smp = iepdf.sample();
-		temp.set_col ( i, smp );
-	}
-
-	return temp;
-}
-
-template<class EPDF>
-double mpdf_internal<EPDF>::evallogcond ( const vec &dt, const vec &cond ) {
-	double tmp;
-	condition ( cond );
-	tmp = iepdf.evallog ( dt );
-	// it_assert_debug(std::isfinite(tmp), "Infinite value");
-	return tmp;
-}
-
-template<class EPDF>
-vec mpdf_internal<EPDF>::evallogcond_m ( const mat &Dt, const vec &cond ) {
-	condition ( cond );
-	return iepdf.evallog_m ( Dt );
-}
-
-template<class EPDF>
-		vec mpdf_internal<EPDF>::evallogcond_m ( const Array<vec> &Dt, const vec &cond ) {
-	condition ( cond );
-	return iepdf.evallog_m ( Dt );
-}
 
 void mpdf::from_setting ( const Setting &set ) {
Index: library/bdm/stat/emix.cpp
===================================================================
--- library/bdm/stat/emix.cpp (revision 477)
+++ library/bdm/stat/emix.cpp (revision 488)
@@ -199,4 +199,19 @@
 	return aprgiw;
 };
+
+vec mmix::samplecond(const vec &cond) {
+	//Sample which component
+	vec cumDist = cumsum ( w );
+	double u0;
+#pragma omp critical
+	u0 = UniRNG.sample();
+
+	int i = 0;
+	while ( ( cumDist ( i ) < u0 ) && ( i < ( w.length() - 1 ) ) ) {
+		i++;
+	}
+
+	return Coms ( i )->samplecond(cond);
+}
 
 }
Index: library/bdm/stat/exp_family.h
===================================================================
--- library/bdm/stat/exp_family.h (revision 487)
+++ library/bdm/stat/exp_family.h (revision 488)
@@ -24,1178 +24,1145 @@
 
 //! Global Uniform_RNG
-	extern Uniform_RNG UniRNG;
+extern Uniform_RNG UniRNG;
 //! Global Normal_RNG
-	extern Normal_RNG NorRNG;
+extern Normal_RNG NorRNG;
 //! Global Gamma_RNG
-	extern Gamma_RNG GamRNG;
-
-	/*!
-	* \brief General conjugate exponential family posterior density.
-
-	* More?...
-	*/
-
-	class eEF : public epdf
-	{
-		public:
+extern Gamma_RNG GamRNG;
+
+/*!
+* \brief General conjugate exponential family posterior density.
+
+* More?...
+*/
+
+class eEF : public epdf
+{
+	public:
 //	eEF() :epdf() {};
-			//! default constructor
-			eEF ( ) :epdf ( ) {};
-			//! logarithm of the normalizing constant, \f$\mathcal{I}\f$
-			virtual double lognc() const =0;
-			//!TODO decide if it is really needed
-			virtual void dupdate ( mat &v ) {it_error ( "Not implemented" );};
-			//!Evaluate normalized log-probability
-			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;}
-			//!Evaluate normalized log-probability for many samples
-			virtual vec evallog_m ( const mat &Val ) const
-			{
-				vec x ( Val.cols() );
-				for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog_nn ( Val.get_col ( i ) ) ;}
-				return x-lognc();
+		//! default constructor
+		eEF () : epdf () {};
+		//! logarithm of the normalizing constant, \f$\mathcal{I}\f$
+		virtual double lognc() const = 0;
+		//!TODO decide if it is really needed
+		virtual void dupdate (mat &v) {it_error ("Not implemented");};
+		//!Evaluate normalized log-probability
+		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;
+		}
+		//!Evaluate normalized log-probability for many samples
+		virtual vec evallog_m (const mat &Val) const {
+			vec x (Val.cols());
+			for (int i = 0;i < Val.cols();i++) {x (i) = evallog_nn (Val.get_col (i)) ;}
+			return x -lognc();
+		}
+		//!Evaluate normalized log-probability for many samples
+		virtual vec evallog_m (const Array<vec> &Val) const {
+			vec x (Val.length());
+			for (int i = 0;i < Val.length();i++) {x (i) = evallog_nn (Val (i)) ;}
+			return x -lognc();
+		}
+		//!Power of the density, used e.g. to flatten the density
+		virtual void pow (double p) {it_error ("Not implemented");};
+};
+
+
+//! Estimator for Exponential family
+class BMEF : public BM
+{
+	protected:
+		//! forgetting factor
+		double frg;
+		//! cached value of lognc() in the previous step (used in evaluation of \c ll )
+		double last_lognc;
+	public:
+		//! Default constructor (=empty constructor)
+		BMEF (double frg0 = 1.0) : BM (), frg (frg0) {}
+		//! Copy constructor
+		BMEF (const BMEF &B) : BM (B), frg (B.frg), last_lognc (B.last_lognc) {}
+		//!get statistics from another model
+		virtual void set_statistics (const BMEF* BM0) {it_error ("Not implemented");};
+		//! Weighted update of sufficient statistics (Bayes rule)
+		virtual void bayes (const vec &data, const double w) {};
+		//original Bayes
+		void bayes (const vec &dt);
+		//!Flatten the posterior according to the given BMEF (of the same type!)
+		virtual void flatten (const BMEF * B) {it_error ("Not implemented");}
+		//!Flatten the posterior as if to keep nu0 data
+//	virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
+
+		BMEF* _copy_ () const {it_error ("function _copy_ not implemented for this BM"); return NULL;};
+};
+
+template<class sq_T, template <typename> class TEpdf >
+class mlnorm;
+
+/*!
+* \brief Gaussian density with positive definite (decomposed) covariance matrix.
+
+* More?...
+*/
+template<class sq_T>
+class enorm : public eEF
+{
+	protected:
+		//! mean value
+		vec mu;
+		//! Covariance matrix in decomposed form
+		sq_T R;
+	public:
+		//!\name Constructors
+		//!@{
+
+		enorm () : eEF (), mu (), R () {};
+		enorm (const vec &mu, const sq_T &R) {set_parameters (mu, R);}
+		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();
+		}
+		//!@}
+
+		//! \name Mathematical operations
+		//!@{
+
+		//! dupdate in exponential form (not really handy)
+		void dupdate (mat &v, double nu = 1.0);
+
+		vec sample() const;
+
+		double evallog_nn (const vec &val) const;
+		double lognc () const;
+		vec mean() const {return mu;}
+		vec variance() const {return diag (R.to_mat());}
+//	mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
+		mpdf* condition (const RV &rvn) const ;
+		enorm<sq_T>* marginal (const RV &rv) const;
+//			epdf* marginal ( const RV &rv ) const;
+		//!@}
+
+		//! \name Access to attributes
+		//!@{
+
+		vec& _mu() {return mu;}
+		void set_mu (const vec mu0) { mu = mu0;}
+		sq_T& _R() {return R;}
+		const sq_T& _R() const {return R;}
+		//!@}
+
+};
+UIREGISTER (enorm<chmat>);
+UIREGISTER (enorm<ldmat>);
+UIREGISTER (enorm<fsqmat>);
+
+
+/*!
+* \brief Gauss-inverse-Wishart density stored in LD form
+
+* For \f$p\f$-variate densities, given rv.count() should be \f$p\times\f$ V.rows().
+*
+*/
+class egiw : public eEF
+{
+	protected:
+		//! Extended information matrix of sufficient statistics
+		ldmat V;
+		//! Number of data records (degrees of freedom) of sufficient statistics
+		double nu;
+		//! Dimension of the output
+		int dimx;
+		//! Dimension of the regressor
+		int nPsi;
+	public:
+		//!\name Constructors
+		//!@{
+		egiw() : eEF() {};
+		egiw (int dimx0, ldmat V0, double nu0 = -1.0) : eEF() {set_parameters (dimx0, V0, nu0);};
+
+		void set_parameters (int dimx0, ldmat V0, double nu0 = -1.0) {
+			dimx = dimx0;
+			nPsi = V0.rows() - dimx;
+			dim = dimx * (dimx + nPsi); // size(R) + size(Theta)
+
+			V = V0;
+			if (nu0 < 0) {
+				nu = 0.1 + nPsi + 2 * dimx + 2; // +2 assures finite expected value of R
+				// terms before that are sufficient for finite normalization
+			} else {
+				nu = nu0;
 			}
-			//!Power of the density, used e.g. to flatten the density
-			virtual void pow ( double p ) {it_error ( "Not implemented" );};
-	};
-
-
-//! Estimator for Exponential family
-	class BMEF : public BM
-	{
-		protected:
-			//! forgetting factor
-			double frg;
-			//! cached value of lognc() in the previous step (used in evaluation of \c ll )
-			double last_lognc;
-		public:
-			//! Default constructor (=empty constructor)
-			BMEF ( double frg0=1.0 ) :BM ( ), frg ( frg0 ) {}
-			//! Copy constructor
-			BMEF ( const BMEF &B ) :BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
-			//!get statistics from another model
-			virtual void set_statistics ( const BMEF* BM0 ) {it_error ( "Not implemented" );};
-			//! Weighted update of sufficient statistics (Bayes rule)
-			virtual void bayes ( const vec &data, const double w ) {};
-			//original Bayes
-			void bayes ( const vec &dt );
-			//!Flatten the posterior according to the given BMEF (of the same type!)
-			virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
-			//!Flatten the posterior as if to keep nu0 data
-//	virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}
-
-			BMEF* _copy_ () const {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
-	};
-
-	template<class sq_T, template <typename> class TEpdf >
-	class mlnorm;
-
-	/*!
-	* \brief Gaussian density with positive definite (decomposed) covariance matrix.
-
-	* More?...
-	*/
-	template<class sq_T>
-	class enorm : public eEF
-	{
-		protected:
-			//! mean value
-			vec mu;
-			//! Covariance matrix in decomposed form
-			sq_T R;
-		public:
-			//!\name Constructors
-			//!@{
-
-			enorm ( ) :eEF ( ), mu ( ),R ( ) {};
-			enorm ( const vec &mu,const sq_T &R ) {set_parameters ( mu,R );}
-			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();
+		}
+		//!@}
+
+		vec sample() const;
+		vec mean() const;
+		vec variance() const;
+
+		//! LS estimate of \f$\theta\f$
+		vec est_theta() const;
+
+		//! Covariance of the LS estimate
+		ldmat est_theta_cov() const;
+
+		void mean_mat (mat &M, mat&R) const;
+		//! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
+		double evallog_nn (const vec &val) const;
+		double lognc () const;
+		void pow (double p) {V *= p;nu *= p;};
+
+		//! \name Access attributes
+		//!@{
+
+		ldmat& _V() {return V;}
+		const ldmat& _V() const {return V;}
+		double& _nu()  {return nu;}
+		const double& _nu() const {return nu;}
+		void from_setting (const Setting &set) {
+			UI::get (nu, set, "nu", UI::compulsory);
+			UI::get (dimx, set, "dimx", UI::compulsory);
+			mat V;
+			UI::get (V, set, "V", UI::compulsory);
+			set_parameters (dimx, V, nu);
+			RV* rv = UI::build<RV> (set, "rv", UI::compulsory);
+			set_rv (*rv);
+			delete rv;
+		}
+		//!@}
+};
+UIREGISTER (egiw);
+
+/*! \brief Dirichlet posterior density
+
+Continuous Dirichlet density of \f$n\f$-dimensional variable \f$x\f$
+\f[
+f(x|\beta) = \frac{\Gamma[\gamma]}{\prod_{i=1}^{n}\Gamma(\beta_i)} \prod_{i=1}^{n}x_i^{\beta_i-1}
+\f]
+where \f$\gamma=\sum_i \beta_i\f$.
+*/
+class eDirich: public eEF
+{
+	protected:
+		//!sufficient statistics
+		vec beta;
+	public:
+		//!\name Constructors
+		//!@{
+
+		eDirich () : eEF () {};
+		eDirich (const eDirich &D0) : eEF () {set_parameters (D0.beta);};
+		eDirich (const vec &beta0) {set_parameters (beta0);};
+		void set_parameters (const vec &beta0) {
+			beta = beta0;
+			dim = beta.length();
+		}
+		//!@}
+
+		vec sample() const {it_error ("Not implemented");return vec_1 (0.0);};
+		vec mean() const {return beta / sum (beta);};
+		vec variance() const {double gamma = sum (beta); return elem_mult (beta, (beta + 1)) / (gamma* (gamma + 1));}
+		//! In this instance, val is ...
+		double evallog_nn (const vec &val) const {
+			double tmp; tmp = (beta - 1) * log (val);
+// 				it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
+			return tmp;
+		};
+		double lognc () const {
+			double tmp;
+			double gam = sum (beta);
+			double lgb = 0.0;
+			for (int i = 0;i < beta.length();i++) {lgb += lgamma (beta (i));}
+			tmp = lgb - lgamma (gam);
+// 				it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
+			return tmp;
+		};
+		//!access function
+		vec& _beta()  {return beta;}
+		//!Set internal parameters
+};
+
+//! \brief Estimator for Multinomial density
+class multiBM : public BMEF
+{
+	protected:
+		//! Conjugate prior and posterior
+		eDirich est;
+		//! Pointer inside est to sufficient statistics
+		vec &beta;
+	public:
+		//!Default constructor
+		multiBM () : BMEF (), est (), beta (est._beta()) {
+			if (beta.length() > 0) {last_lognc = est.lognc();}
+			else{last_lognc = 0.0;}
+		}
+		//!Copy constructor
+		multiBM (const multiBM &B) : BMEF (B), est (B.est), beta (est._beta()) {}
+		//! Sets sufficient statistics to match that of givefrom mB0
+		void set_statistics (const BM* mB0) {const multiBM* mB = dynamic_cast<const multiBM*> (mB0); beta = mB->beta;}
+		void bayes (const vec &dt) {
+			if (frg < 1.0) {beta *= frg;last_lognc = est.lognc();}
+			beta += dt;
+			if (evalll) {ll = est.lognc() - last_lognc;}
+		}
+		double logpred (const vec &dt) const {
+			eDirich pred (est);
+			vec &beta = pred._beta();
+
+			double lll;
+			if (frg < 1.0)
+				{beta *= frg;lll = pred.lognc();}
+			else
+				if (evalll) {lll = last_lognc;}
+				else{lll = pred.lognc();}
+
+			beta += dt;
+			return pred.lognc() - lll;
+		}
+		void flatten (const BMEF* B) {
+			const multiBM* E = dynamic_cast<const multiBM*> (B);
+			// sum(beta) should be equal to sum(B.beta)
+			const vec &Eb = E->beta;//const_cast<multiBM*> ( E )->_beta();
+			beta *= (sum (Eb) / sum (beta));
+			if (evalll) {last_lognc = est.lognc();}
+		}
+		const epdf& posterior() const {return est;};
+		const eDirich* _e() const {return &est;};
+		void set_parameters (const vec &beta0) {
+			est.set_parameters (beta0);
+			if (evalll) {last_lognc = est.lognc();}
+		}
+};
+
+/*!
+ \brief Gamma posterior density
+
+ Multivariate Gamma density as product of independent univariate densities.
+ \f[
+ f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i)
+ \f]
+*/
+
+class egamma : public eEF
+{
+	protected:
+		//! Vector \f$\alpha\f$
+		vec alpha;
+		//! Vector \f$\beta\f$
+		vec beta;
+	public :
+		//! \name Constructors
+		//!@{
+		egamma () : eEF (), alpha (0), beta (0) {};
+		egamma (const vec &a, const vec &b) {set_parameters (a, b);};
+		void set_parameters (const vec &a, const vec &b) {alpha = a, beta = b;dim = alpha.length();};
+		//!@}
+
+		vec sample() const;
+		//! TODO: is it used anywhere?
+//	mat sample ( int N ) const;
+		double evallog (const vec &val) const;
+		double lognc () const;
+		//! Returns poiter to alpha and beta. Potentially dengerous: use with care!
+		vec& _alpha() {return alpha;}
+		vec& _beta() {return beta;}
+		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::compulsory);
+			UI::get (beta, set, "beta", UI::compulsory);
+			validate();
+		}
+		void validate() {
+			it_assert (alpha.length() == beta.length(), "parameters do not match");
+			dim = alpha.length();
+		}
+};
+UIREGISTER (egamma);
+/*!
+ \brief Inverse-Gamma posterior density
+
+ Multivariate inverse-Gamma density as product of independent univariate densities.
+ \f[
+ f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i)
+ \f]
+
+Vector \f$\beta\f$ has different meaning (in fact it is 1/beta as used in definition of iG)
+
+ Inverse Gamma can be converted to Gamma using \f[
+ x\sim iG(a,b) => 1/x\sim G(a,1/b)
+\f]
+This relation is used in sampling.
+ */
+
+class eigamma : public egamma
+{
+	protected:
+	public :
+		//! \name Constructors
+		//! All constructors are inherited
+		//!@{
+		//!@}
+
+		vec sample() const {return 1.0 / egamma::sample();};
+		//! Returns poiter to alpha and beta. Potentially dangerous: use with care!
+		vec mean() const {return elem_div (beta, alpha - 1);}
+		vec variance() const {vec mea = mean(); return elem_div (elem_mult (mea, mea), alpha - 2);}
+};
+/*
+//! Weighted mixture of epdfs with external owned components.
+class emix : public epdf {
+protected:
+	int n;
+	vec &w;
+	Array<epdf*> Coms;
+public:
+//! Default constructor
+	emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
+	void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
+	vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
+	vec sample() {it_error ( "Not implemented" );return 0;}
+};
+*/
+
+//!  Uniform distributed density on a rectangular support
+
+class euni: public epdf
+{
+	protected:
+//! lower bound on support
+		vec low;
+//! upper bound on support
+		vec high;
+//! internal
+		vec distance;
+//! normalizing coefficients
+		double nk;
+//! cache of log( \c nk )
+		double lnk;
+	public:
+		//! \name Constructors
+		//!@{
+		euni () : epdf () {}
+		euni (const vec &low0, const vec &high0) {set_parameters (low0, high0);}
+		void set_parameters (const vec &low0, const vec &high0) {
+			distance = high0 - low0;
+			it_assert_debug (min (distance) > 0.0, "bad support");
+			low = low0;
+			high = high0;
+			nk = prod (1.0 / distance);
+			lnk = log (nk);
+			dim = low.length();
+		}
+		//!@}
+
+		double eval (const vec &val) const  {return nk;}
+		double evallog (const vec &val) const  {
+			if (any (val < low) && any (val > high)) {return inf;}
+			else return lnk;
+		}
+		vec sample() const {
+			vec smp (dim);
+#pragma omp critical
+			UniRNG.sample_vector (dim , smp);
+			return low + elem_mult (distance, smp);
+		}
+		//! set values of \c low and \c high
+		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::compulsory);
+			UI::get (low, set, "low", UI::compulsory);
+		}
+};
+
+
+/*!
+ \brief Normal distributed linear function with linear function of mean value;
+
+ Mean value \f$mu=A*rvc+mu_0\f$.
+*/
+template < class sq_T, template <typename> class TEpdf = enorm >
+class mlnorm : public mpdf_internal< TEpdf<sq_T> >
+{
+	protected:
+		//! Internal epdf that arise by conditioning on \c rvc
+		mat A;
+		vec mu_const;
+//			vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
+	public:
+		//! \name Constructors
+		//!@{
+		mlnorm() : mpdf_internal< TEpdf<sq_T> >() {};
+		mlnorm (const mat &A, const vec &mu0, const sq_T &R) : mpdf_internal< TEpdf<sq_T> >() {
+			set_parameters (A, mu0, R);
+		}
+
+		//! Set \c A and \c R
+		void set_parameters (const  mat &A0, const vec &mu0, const sq_T &R0) {
+			it_assert_debug (A0.rows() == mu0.length(), "");
+			it_assert_debug (A0.rows() == R0.rows(), "");
+
+			this->iepdf.set_parameters (zeros (A0.rows()), R0);
+			A = A0;
+			mu_const = mu0;
+			this->dimc = A0.cols();
+		}
+		//!@}
+		//! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
+		void condition (const vec &cond) {
+			this->iepdf._mu() = A * cond + mu_const;
+//R is already assigned;
+		}
+
+		//!access function
+		vec& _mu_const() {return mu_const;}
+		//!access function
+		mat& _A() {return A;}
+		//!access function
+		mat _R() { return this->iepdf._R().to_mat(); }
+
+		template<typename sq_M>
+		friend std::ostream &operator<< (std::ostream &os,  mlnorm<sq_M, enorm> &ml);
+
+		void from_setting (const Setting &set) {
+			mpdf::from_setting (set);
+
+			UI::get (A, set, "A", UI::compulsory);
+			UI::get (mu_const, set, "const", UI::compulsory);
+			mat R0;
+			UI::get (R0, set, "R", UI::compulsory);
+			set_parameters (A, mu_const, R0);
+		};
+};
+UIREGISTER (mlnorm<ldmat>);
+UIREGISTER (mlnorm<fsqmat>);
+UIREGISTER (mlnorm<chmat>);
+
+//! Mpdf with general function for mean value
+template<class sq_T>
+class mgnorm : public mpdf_internal< enorm< sq_T > >
+{
+	protected:
+//			vec &mu; WHY NOT?
+		fnc* g;
+	public:
+		//!default constructor
+		mgnorm() : mpdf_internal<enorm<sq_T> >() { }
+		//!set mean function
+		inline void set_parameters (fnc* g0, const sq_T &R0);
+		inline void condition (const vec &cond);
+
+
+		/*! UI for mgnorm
+
+		The mgnorm is constructed from a structure with fields:
+		\code
+		system = {
+			type = "mgnorm";
+			// function for mean value evolution
+			g = {type="fnc"; ... }
+
+			// variance
+			R = [1, 0,
+				 0, 1];
+			// --OR --
+			dR = [1, 1];
+
+			// == OPTIONAL ==
+
+			// description of y variables
+			y = {type="rv"; names=["y", "u"];};
+			// description of u variable
+			u = {type="rv"; names=[];}
+		};
+		\endcode
+
+		Result if
+		*/
+
+		void from_setting (const Setting &set) {
+			fnc* g = UI::build<fnc> (set, "g", UI::compulsory);
+
+			mat R;
+			vec dR;
+			if (UI::get (dR, set, "dR"))
+				R = diag (dR);
+			else
+				UI::get (R, set, "R", UI::compulsory);
+
+			set_parameters (g, R);
+		}
+};
+
+UIREGISTER (mgnorm<chmat>);
+
+
+/*! (Approximate) Student t density with linear function of mean value
+
+The internal epdf of this class is of the type of a Gaussian (enorm).
+However, each conditioning is trying to assure the best possible approximation by taking into account the zeta function. See [] for reference.
+
+Perhaps a moment-matching technique?
+*/
+class mlstudent : public mlnorm<ldmat, enorm>
+{
+	protected:
+		ldmat Lambda;
+		ldmat &_R;
+		ldmat Re;
+	public:
+		mlstudent () : mlnorm<ldmat, enorm> (),
+				Lambda (),	_R (iepdf._R()) {}
+		void set_parameters (const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0) {
+			it_assert_debug (A0.rows() == mu0.length(), "");
+			it_assert_debug (R0.rows() == A0.rows(), "");
+
+			iepdf.set_parameters (mu0, Lambda);  //
+			A = A0;
+			mu_const = mu0;
+			Re = R0;
+			Lambda = Lambda0;
+		}
+		void condition (const vec &cond) {
+			iepdf._mu() = A * cond + mu_const;
+			double zeta;
+			//ugly hack!
+			if ( (cond.length() + 1) == Lambda.rows()) {
+				zeta = Lambda.invqform (concat (cond, vec_1 (1.0)));
+			} else {
+				zeta = Lambda.invqform (cond);
 			}
-			//!@}
-
-			//! \name Mathematical operations
-			//!@{
-
-			//! dupdate in exponential form (not really handy)
-			void dupdate ( mat &v,double nu=1.0 );
-
-			vec sample() const;
-
-			double evallog_nn ( const vec &val ) const;
-			double lognc () const;
-			vec mean() const {return mu;}
-			vec variance() const {return diag ( R.to_mat() );}
-//	mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
-			mpdf* condition ( const RV &rvn ) const ;
-			enorm<sq_T>* marginal ( const RV &rv ) const;
-//			epdf* marginal ( const RV &rv ) const;
-			//!@}
-
-			//! \name Access to attributes
-			//!@{
-
-			vec& _mu() {return mu;}
-			void set_mu ( const vec mu0 ) { mu=mu0;}
-			sq_T& _R() {return R;}
-			const sq_T& _R() const {return R;}
-			//!@}
-
-	};
-	UIREGISTER(enorm<chmat>);
-	UIREGISTER(enorm<ldmat>);
-	UIREGISTER(enorm<fsqmat>);
-
-
-	/*!
-	* \brief Gauss-inverse-Wishart density stored in LD form
-
-	* For \f$p\f$-variate densities, given rv.count() should be \f$p\times\f$ V.rows().
-	*
-	*/
-	class egiw : public eEF
-	{
-		protected:
-			//! Extended information matrix of sufficient statistics
-			ldmat V;
-			//! Number of data records (degrees of freedom) of sufficient statistics
-			double nu;
-			//! Dimension of the output
-			int dimx;
-			//! Dimension of the regressor
-			int nPsi;
-		public:
-			//!\name Constructors
-			//!@{
-			egiw() :eEF() {};
-			egiw ( int dimx0, ldmat V0, double nu0=-1.0 ) :eEF() {set_parameters ( dimx0,V0, nu0 );};
-
-			void set_parameters ( int dimx0, ldmat V0, double nu0=-1.0 )
-			{
-				dimx=dimx0;
-				nPsi = V0.rows()-dimx;
-				dim = dimx* ( dimx+nPsi ); // size(R) + size(Theta)
-
-				V=V0;
-				if ( nu0<0 )
-				{
-					nu = 0.1 +nPsi +2*dimx +2; // +2 assures finite expected value of R
-					// terms before that are sufficient for finite normalization
-				}
-				else
-				{
-					nu=nu0;
+			_R = Re;
+			_R *= (1 + zeta);// / ( nu ); << nu is in Re!!!!!!
+		};
+
+};
+/*!
+\brief  Gamma random walk
+
+Mean value, \f$\mu\f$, of this density is given by \c rvc .
+Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
+This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
+
+The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
+*/
+class mgamma : public mpdf_internal<egamma>
+{
+	protected:
+
+		//! Constant \f$k\f$
+		double k;
+
+		//! cache of iepdf.beta
+		vec &_beta;
+
+	public:
+		//! Constructor
+		mgamma() : mpdf_internal<egamma>(), k (0),
+				_beta (iepdf._beta()) {
+		}
+
+		//! Set value of \c k
+		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", UI::compulsory);
+			UI::get (k, set, "k", UI::compulsory);
+			set_parameters (k, betatmp);
+		}
+};
+UIREGISTER (mgamma);
+
+/*!
+\brief  Inverse-Gamma random walk
+
+Mean value, \f$ \mu \f$, of this density is given by \c rvc .
+Standard deviation of the random walk is proportional to one \f$ k \f$-th the mean.
+This is achieved by setting \f$ \alpha=\mu/k^2+2 \f$ and \f$ \beta=\mu(\alpha-1)\f$.
+
+The standard deviation of the walk is then: \f$ \mu/\sqrt(k)\f$.
+ */
+class migamma : public mpdf_internal<eigamma>
+{
+	protected:
+		//! Constant \f$k\f$
+		double k;
+
+		//! cache of iepdf.alpha
+		vec &_alpha;
+
+		//! cache of iepdf.beta
+		vec &_beta;
+
+	public:
+		//! \name Constructors
+		//!@{
+		migamma() : mpdf_internal<eigamma>(),
+				k (0),
+				_alpha (iepdf._alpha()),
+				_beta (iepdf._beta()) {
+		}
+
+		migamma (const migamma &m) : mpdf_internal<eigamma>(),
+				k (0),
+				_alpha (iepdf._alpha()),
+				_beta (iepdf._beta()) {
+		}
+		//!@}
+
+		//! Set value of \c k
+		void set_parameters (int len, double k0) {
+			k = k0;
+			iepdf.set_parameters ( (1.0 / (k*k) + 2.0) *ones (len) /*alpha*/, ones (len) /*beta*/);
+			dimc = dimension();
+		};
+		void condition (const vec &val) {
+			_beta = elem_mult (val, (_alpha - 1.0));
+		};
+};
+
+
+/*!
+\brief  Gamma random walk around a fixed point
+
+Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
+\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
+
+Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
+This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
+
+The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
+*/
+class mgamma_fix : public mgamma
+{
+	protected:
+		//! parameter l
+		double l;
+		//! reference vector
+		vec refl;
+	public:
+		//! Constructor
+		mgamma_fix () : mgamma (), refl () {};
+		//! Set value of \c k
+		void set_parameters (double k0 , vec ref0, double l0) {
+			mgamma::set_parameters (k0, ref0);
+			refl = pow (ref0, 1.0 - l0);l = l0;
+			dimc = dimension();
+		};
+
+		void condition (const vec &val) {vec mean = elem_mult (refl, pow (val, l)); _beta = k / mean;};
+};
+
+
+/*!
+\brief  Inverse-Gamma random walk around a fixed point
+
+Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
+\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
+
+==== Check == vv =
+Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
+This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
+
+The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
+ */
+class migamma_ref : public migamma
+{
+	protected:
+		//! parameter l
+		double l;
+		//! reference vector
+		vec refl;
+	public:
+		//! Constructor
+		migamma_ref () : migamma (), refl () {};
+		//! Set value of \c k
+		void set_parameters (double k0 , vec ref0, double l0) {
+			migamma::set_parameters (ref0.length(), k0);
+			refl = pow (ref0, 1.0 - l0);
+			l = l0;
+			dimc = dimension();
+		};
+
+		void condition (const vec &val) {
+			vec mean = elem_mult (refl, pow (val, l));
+			migamma::condition (mean);
+		};
+
+		/*! UI for migamma_ref
+
+		The migamma_ref is constructed from a structure with fields:
+		\code
+		system = {
+			type = "migamma_ref";
+			ref = [1e-5; 1e-5; 1e-2 1e-3];            // reference vector
+			l = 0.999;                                // constant l
+			k = 0.1;                                  // constant k
+			
+			// == OPTIONAL ==
+			// description of y variables
+			y = {type="rv"; names=["y", "u"];};
+			// description of u variable
+			u = {type="rv"; names=[];}
+		};
+		\endcode
+
+		Result if
+		 */
+		void from_setting (const Setting &set);
+
+		// TODO dodelat void to_setting( Setting &set ) const;
+};
+
+
+UIREGISTER (migamma_ref);
+
+/*! Log-Normal probability density
+ only allow diagonal covariances!
+
+Density of the form \f$ \log(x)\sim \mathcal{N}(\mu,\sigma^2) \f$ , i.e.
+\f[
+x \sim \frac{1}{x\sigma\sqrt{2\pi}}\exp{-\frac{1}{2\sigma^2}(\log(x)-\mu)}
+\f]
+
+*/
+class elognorm: public enorm<ldmat>
+{
+	public:
+		vec sample() const {return exp (enorm<ldmat>::sample());};
+		vec mean() const {vec var = enorm<ldmat>::variance();return exp (mu - 0.5*var);};
+
+};
+
+/*!
+\brief  Log-Normal random walk
+
+Mean value, \f$\mu\f$, is...
+
+==== Check == vv =
+Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
+This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
+
+The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
+ */
+class mlognorm : public mpdf_internal<elognorm>
+{
+	protected:
+		//! parameter 1/2*sigma^2
+		double sig2;
+
+		//! access
+		vec &mu;
+	public:
+		//! Constructor
+		mlognorm() : mpdf_internal<elognorm>(),
+				sig2 (0),
+				mu (iepdf._mu()) {
+		}
+
+		//! Set value of \c k
+		void set_parameters (int size, double k) {
+			sig2 = 0.5 * log (k * k + 1);
+			iepdf.set_parameters (zeros (size), 2*sig2*eye (size));
+
+			dimc = size;
+		};
+
+		void condition (const vec &val) {
+			mu = log (val) - sig2;//elem_mult ( refl,pow ( val,l ) );
+		};
+
+		/*! UI for mlognorm
+
+		The mlognorm is constructed from a structure with fields:
+		\code
+		system = {
+			type = "mlognorm";
+			k = 0.1;                                  // constant k
+			mu0 = [1., 1.];
+			
+			// == OPTIONAL ==
+			// description of y variables
+			y = {type="rv"; names=["y", "u"];};
+			// description of u variable
+			u = {type="rv"; names=[];}
+		};
+		\endcode
+
+		 */
+		void from_setting (const Setting &set);
+
+		// TODO dodelat void to_setting( Setting &set ) const;
+
+};
+
+UIREGISTER (mlognorm);
+
+/*! inverse Wishart density defined on Choleski decomposition
+
+*/
+class eWishartCh : public epdf
+{
+	protected:
+		//! Upper-Triagle of Choleski decomposition of \f$ \Psi \f$
+		chmat Y;
+		//! dimension of matrix  \f$ \Psi \f$
+		int p;
+		//! degrees of freedom  \f$ \nu \f$
+		double delta;
+	public:
+		void set_parameters (const mat &Y0, const double delta0) {Y = chmat (Y0);delta = delta0; p = Y.rows(); dim = p * p; }
+		mat sample_mat() const {
+			mat X = zeros (p, p);
+
+			//sample diagonal
+			for (int i = 0;i < p;i++) {
+				GamRNG.setup (0.5* (delta - i) , 0.5);   // no +1 !! index if from 0
+#pragma omp critical
+				X (i, i) = sqrt (GamRNG());
+			}
+			//do the rest
+			for (int i = 0;i < p;i++) {
+				for (int j = i + 1;j < p;j++) {
+#pragma omp critical
+					X (i, j) = NorRNG.sample();
 				}
 			}
-			//!@}
-
-			vec sample() const;
-			vec mean() const;
-			vec variance() const;
-
-			//! LS estimate of \f$\theta\f$
-			vec est_theta() const;
-
-			//! Covariance of the LS estimate
-			ldmat est_theta_cov() const;
-
-			void mean_mat ( mat &M, mat&R ) const;
-			//! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
-			double evallog_nn ( const vec &val ) const;
-			double lognc () const;
-			void pow ( double p ) {V*=p;nu*=p;};
-
-			//! \name Access attributes
-			//!@{
-
-			ldmat& _V() {return V;}
-			const ldmat& _V() const {return V;}
-			double& _nu()  {return nu;}
-			const double& _nu() const {return nu;}
-			void from_setting(const Setting &set){
-				UI::get(nu, set, "nu", UI::compulsory);
-				UI::get(dimx, set, "dimx", UI::compulsory);
-				mat V;
-				UI::get(V,set,"V", UI::compulsory);
-				set_parameters(dimx, V, nu);
-				RV* rv=UI::build<RV>(set,"rv", UI::compulsory);
-				set_rv(*rv);
-				delete rv;
+			return X*Y._Ch();// return upper triangular part of the decomposition
+		}
+		vec sample () const {
+			return vec (sample_mat()._data(), p*p);
+		}
+		//! fast access function y0 will be copied into Y.Ch.
+		void setY (const mat &Ch0) {copy_vector (dim, Ch0._data(), Y._Ch()._data());}
+		//! fast access function y0 will be copied into Y.Ch.
+		void _setY (const vec &ch0) {copy_vector (dim, ch0._data(), Y._Ch()._data()); }
+		//! access function
+		const chmat& getY() const {return Y;}
+};
+
+class eiWishartCh: public epdf
+{
+	protected:
+		eWishartCh W;
+		int p;
+		double delta;
+	public:
+		void set_parameters (const mat &Y0, const double delta0) {
+			delta = delta0;
+			W.set_parameters (inv (Y0), delta0);
+			dim = W.dimension(); p = Y0.rows();
+		}
+		vec sample() const {mat iCh; iCh = inv (W.sample_mat()); return vec (iCh._data(), dim);}
+		void _setY (const vec &y0) {
+			mat Ch (p, p);
+			mat iCh (p, p);
+			copy_vector (dim, y0._data(), Ch._data());
+
+			iCh = inv (Ch);
+			W.setY (iCh);
+		}
+		virtual double evallog (const vec &val) const {
+			chmat X (p);
+			const chmat& Y = W.getY();
+
+			copy_vector (p*p, val._data(), X._Ch()._data());
+			chmat iX (p);X.inv (iX);
+			// compute
+// 				\frac{ |\Psi|^{m/2}|X|^{-(m+p+1)/2}e^{-tr(\Psi X^{-1})/2} }{ 2^{mp/2}\Gamma_p(m/2)},
+			mat M = Y.to_mat() * iX.to_mat();
+
+			double log1 = 0.5 * p * (2 * Y.logdet()) - 0.5 * (delta + p + 1) * (2 * X.logdet()) - 0.5 * trace (M);
+			//Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!
+
+			/*				if (0) {
+								mat XX=X.to_mat();
+								mat YY=Y.to_mat();
+
+								double log2 = 0.5*p*log(det(YY))-0.5*(delta+p+1)*log(det(XX))-0.5*trace(YY*inv(XX));
+								cout << log1 << "," << log2 << endl;
+							}*/
+			return log1;
+		};
+
+};
+
+class rwiWishartCh : public mpdf
+{
+	protected:
+		eiWishartCh eiW;
+		//!square root of \f$ \nu-p-1 \f$ - needed for computation of \f$ \Psi \f$ from conditions
+		double sqd;
+		//reference point for diagonal
+		vec refl;
+		double l;
+		int p;
+
+	public:
+		rwiWishartCh() : eiW(),
+				sqd (0), l (0), p (0) {
+			set_ep (eiW);
+		}
+
+		void set_parameters (int p0, double k, vec ref0, double l0) {
+			p = p0;
+			double delta = 2 / (k * k) + p + 3;
+			sqd = sqrt (delta - p - 1);
+			l = l0;
+			refl = pow (ref0, 1 - l);
+
+			eiW.set_parameters (eye (p), delta);
+			dimc = eiW.dimension();
+		}
+		void condition (const vec &c) {
+			vec z = c;
+			int ri = 0;
+			for (int i = 0;i < p*p;i += (p + 1)) {//trace diagonal element
+				z (i) = pow (z (i), l) * refl (ri);
+				ri++;
 			}
-			//!@}
-	};
-	UIREGISTER(egiw);
-
-	/*! \brief Dirichlet posterior density
-
-	Continuous Dirichlet density of \f$n\f$-dimensional variable \f$x\f$
-	\f[
-	f(x|\beta) = \frac{\Gamma[\gamma]}{\prod_{i=1}^{n}\Gamma(\beta_i)} \prod_{i=1}^{n}x_i^{\beta_i-1}
-	\f]
-	where \f$\gamma=\sum_i \beta_i\f$.
-	*/
-	class eDirich: public eEF
-	{
-		protected:
-			//!sufficient statistics
-			vec beta;
-		public:
-			//!\name Constructors
-			//!@{
-
-			eDirich () : eEF ( ) {};
-			eDirich ( const eDirich &D0 ) : eEF () {set_parameters ( D0.beta );};
-			eDirich ( const vec &beta0 ) {set_parameters ( beta0 );};
-			void set_parameters ( const vec &beta0 )
-			{
-				beta= beta0;
-				dim = beta.length();
-			}
-			//!@}
-
-			vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
-			vec mean() const {return beta/sum(beta);};
-			vec variance() const {double gamma =sum(beta); return elem_mult ( beta, ( beta+1 ) ) / ( gamma* ( gamma+1 ) );}
-			//! In this instance, val is ...
-			double evallog_nn ( const vec &val ) const
-			{
-				double tmp; tmp= ( beta-1 ) *log ( val );
-// 				it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
-				return tmp;
-			};
-			double lognc () const
-			{
-				double tmp;
-				double gam=sum ( beta );
-				double lgb=0.0;
-				for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
-				tmp= lgb-lgamma ( gam );
-// 				it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
-				return tmp;
-			};
-			//!access function
-			vec& _beta()  {return beta;}
-			//!Set internal parameters
-	};
-
-//! \brief Estimator for Multinomial density
-	class multiBM : public BMEF
-	{
-		protected:
-			//! Conjugate prior and posterior
-			eDirich est;
-			//! Pointer inside est to sufficient statistics
-			vec &beta;
-		public:
-			//!Default constructor
-			multiBM ( ) : BMEF ( ),est ( ),beta ( est._beta() )
-			{
-				if ( beta.length() >0 ) {last_lognc=est.lognc();}
-				else{last_lognc=0.0;}
-			}
-			//!Copy constructor
-			multiBM ( const multiBM &B ) : BMEF ( B ),est ( B.est ),beta ( est._beta() ) {}
-			//! Sets sufficient statistics to match that of givefrom mB0
-			void set_statistics ( const BM* mB0 ) {const multiBM* mB=dynamic_cast<const multiBM*> ( mB0 ); beta=mB->beta;}
-			void bayes ( const vec &dt )
-			{
-				if ( frg<1.0 ) {beta*=frg;last_lognc=est.lognc();}
-				beta+=dt;
-				if ( evalll ) {ll=est.lognc()-last_lognc;}
-			}
-			double logpred ( const vec &dt ) const
-			{
-				eDirich pred ( est );
-				vec &beta = pred._beta();
-
-				double lll;
-				if ( frg<1.0 )
-					{beta*=frg;lll=pred.lognc();}
-				else
-					if ( evalll ) {lll=last_lognc;}
-					else{lll=pred.lognc();}
-
-				beta+=dt;
-				return pred.lognc()-lll;
-			}
-			void flatten ( const BMEF* B )
-			{
-				const multiBM* E=dynamic_cast<const multiBM*> ( B );
-				// sum(beta) should be equal to sum(B.beta)
-				const vec &Eb=E->beta;//const_cast<multiBM*> ( E )->_beta();
-				beta*= ( sum ( Eb ) /sum ( beta ) );
-				if ( evalll ) {last_lognc=est.lognc();}
-			}
-			const epdf& posterior() const {return est;};
-			const eDirich* _e() const {return &est;};
-			void set_parameters ( const vec &beta0 )
-			{
-				est.set_parameters ( beta0 );
-				if ( evalll ) {last_lognc=est.lognc();}
-			}
-	};
-
-	/*!
-	 \brief Gamma posterior density
-
-	 Multivariate Gamma density as product of independent univariate densities.
-	 \f[
-	 f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i)
-	 \f]
-	*/
-
-	class egamma : public eEF
-	{
-		protected:
-			//! Vector \f$\alpha\f$
-			vec alpha;
-			//! Vector \f$\beta\f$
-			vec beta;
-		public :
-			//! \name Constructors
-			//!@{
-			egamma ( ) :eEF ( ), alpha ( 0 ), beta ( 0 ) {};
-			egamma ( const vec &a, const vec &b ) {set_parameters ( a, b );};
-			void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;dim = alpha.length();};
-			//!@}
-
-			vec sample() const;
-			//! TODO: is it used anywhere?
-//	mat sample ( int N ) const;
-			double evallog ( const vec &val ) const;
-			double lognc () const;
-			//! Returns poiter to alpha and beta. Potentially dengerous: use with care!
-			vec& _alpha() {return alpha;}
-			vec& _beta() {return beta;}
-			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::compulsory);
-				UI::get(beta,set,"beta", UI::compulsory);
-				validate();
-			}
-			void validate(){
-				it_assert(alpha.length() ==beta.length(), "parameters do not match");
-				dim =alpha.length();
-			}
-	};
-UIREGISTER(egamma);
-	/*!
-	 \brief Inverse-Gamma posterior density
-
-	 Multivariate inverse-Gamma density as product of independent univariate densities.
-	 \f[
-	 f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i)
-	 \f]
-
-	Vector \f$\beta\f$ has different meaning (in fact it is 1/beta as used in definition of iG)
-
-	 Inverse Gamma can be converted to Gamma using \f[
-	 x\sim iG(a,b) => 1/x\sim G(a,1/b)
-	\f]
-	This relation is used in sampling.
-	 */
-
-	class eigamma : public egamma
-	{
-		protected:
-		public :
-			//! \name Constructors
-			//! All constructors are inherited
-			//!@{
-			//!@}
-
-			vec sample() const {return 1.0/egamma::sample();};
-			//! Returns poiter to alpha and beta. Potentially dangerous: use with care!
-			vec mean() const {return elem_div ( beta,alpha-1 );}
-			vec variance() const {vec mea=mean(); return elem_div ( elem_mult ( mea,mea ),alpha-2 );}
-	};
-	/*
-	//! Weighted mixture of epdfs with external owned components.
-	class emix : public epdf {
-	protected:
+
+			eiW._setY (sqd*z);
+		}
+};
+
+//! Switch between various resampling methods.
+enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
+/*!
+\brief Weighted empirical density
+
+Used e.g. in particle filters.
+*/
+class eEmp: public epdf
+{
+	protected :
+		//! Number of particles
 		int n;
-		vec &w;
-		Array<epdf*> Coms;
-	public:
-	//! Default constructor
-		emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
-		void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
-		vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
-		vec sample() {it_error ( "Not implemented" );return 0;}
-	};
-	*/
-
-//!  Uniform distributed density on a rectangular support
-
-	class euni: public epdf
-	{
-		protected:
-//! lower bound on support
-			vec low;
-//! upper bound on support
-			vec high;
-//! internal
-			vec distance;
-//! normalizing coefficients
-			double nk;
-//! cache of log( \c nk )
-			double lnk;
-		public:
-			//! \name Constructors
-			//!@{
-			euni ( ) :epdf ( ) {}
-			euni ( const vec &low0, const vec &high0 ) {set_parameters ( low0,high0 );}
-			void set_parameters ( const vec &low0, const vec &high0 )
-			{
-				distance = high0-low0;
-				it_assert_debug ( min ( distance ) >0.0,"bad support" );
-				low = low0;
-				high = high0;
-				nk = prod ( 1.0/distance );
-				lnk = log ( nk );
-				dim = low.length();
-			}
-			//!@}
-
-			double eval ( const vec &val ) const  {return nk;}
-			double evallog ( const vec &val ) const  {
-				if (any(val<low) && any(val>high)) {return inf;}
-				else return lnk;
-			}
-			vec sample() const
-			{
-				vec smp ( dim );
-#pragma omp critical
-				UniRNG.sample_vector ( dim ,smp );
-				return low+elem_mult ( distance,smp );
-			}
-			//! set values of \c low and \c high
-			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::compulsory);
-				UI::get(low,set,"low", UI::compulsory);
-			}
-	};
-
-
-	/*!
-	 \brief Normal distributed linear function with linear function of mean value;
-
-	 Mean value \f$mu=A*rvc+mu_0\f$.
-	*/
-	template<class sq_T, template <typename> class TEpdf =enorm>
-	class mlnorm : public mpdf_internal< TEpdf<sq_T> >
-	{
-		protected:
-			//! Internal epdf that arise by conditioning on \c rvc
-			mat A;
-			vec mu_const;
-//			vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
-		public:
-			//! \name Constructors
-			//!@{
-			mlnorm():mpdf_internal< TEpdf<sq_T> >() {};
-			mlnorm (const mat &A, const vec &mu0, const sq_T &R ) : mpdf_internal< TEpdf<sq_T> >()
-			{
-				set_parameters ( A,mu0,R );
-			}
-
-			//! Set \c A and \c R
-			void set_parameters ( const  mat &A0, const vec &mu0, const sq_T &R0 ){
-				it_assert_debug ( A0.rows() ==mu0.length(),"" );
-				it_assert_debug ( A0.rows() ==R0.rows(),"" );
-
-				this->iepdf.set_parameters(zeros(A0.rows()), R0);
-				A = A0;
-				mu_const = mu0;
-				this->dimc = A0.cols();
-			}
-			//!@}
-			//! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
-			void condition ( const vec &cond ){
-				this->iepdf._mu()= A*cond + mu_const;
-//R is already assigned;
-			}
-
-			//!access function
-			vec& _mu_const() {return mu_const;}
-			//!access function
-			mat& _A() {return A;}
-			//!access function
-			mat _R(){ return this->iepdf._R().to_mat(); }
-		
-			template<typename sq_M>
-			friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M,enorm> &ml );
-			
-			void from_setting(const Setting &set){
-				mpdf::from_setting(set);	
-
-				UI::get(A,set,"A", UI::compulsory);
-				UI::get(mu_const,set,"const", UI::compulsory);
-				mat R0;
-				UI::get(R0,set,"R", UI::compulsory);
-				set_parameters(A,mu_const,R0);
-			};
-	};
-	UIREGISTER(mlnorm<ldmat>);
-	UIREGISTER(mlnorm<fsqmat>);
-	UIREGISTER(mlnorm<chmat>);
-
-//! Mpdf with general function for mean value
-	template<class sq_T>
-	class mgnorm : public mpdf_internal< enorm< sq_T > >
-	{
-		protected:
-//			vec &mu; WHY NOT?
-			fnc* g;
-		public:
-			//!default constructor
-			mgnorm():mpdf_internal<enorm<sq_T> >() { }
-			//!set mean function
-			inline void set_parameters ( fnc* g0, const sq_T &R0 );
-			inline void condition ( const vec &cond );
-
-
-			/*! UI for mgnorm
-
-			The mgnorm is constructed from a structure with fields:
-			\code
-			system = {
-				type = "mgnorm";
-				// function for mean value evolution
-				g = {type="fnc"; ... }
-
-				// variance
-				R = [1, 0,
-					 0, 1];
-				// --OR --
-				dR = [1, 1];
-
-				// == OPTIONAL ==
-
-				// description of y variables
-				y = {type="rv"; names=["y", "u"];};
-				// description of u variable
-				u = {type="rv"; names=[];}
-			};
-			\endcode
-
-			Result if
-			*/
-
-			void from_setting( const Setting &set ) 
-			{	
-				fnc* g = UI::build<fnc>( set, "g", UI::compulsory );
-
-				mat R;					
-				vec dR;					
-				if ( UI::get( dR, set, "dR" ) )
-					R=diag(dR);
-				else 
-					UI::get( R, set, "R", UI::compulsory); 
-		
-				set_parameters(g,R);
-			}
-	};
-
-	UIREGISTER(mgnorm<chmat>);
-
-
-	/*! (Approximate) Student t density with linear function of mean value
-
-	The internal epdf of this class is of the type of a Gaussian (enorm).
-	However, each conditioning is trying to assure the best possible approximation by taking into account the zeta function. See [] for reference.
-
-	Perhaps a moment-matching technique?
-	*/
-	class mlstudent : public mlnorm<ldmat,enorm>
-	{
-		protected:
-			ldmat Lambda;
-			ldmat &_R;
-			ldmat Re;
-		public:
-			mlstudent ( ) :mlnorm<ldmat,enorm> (),
-					Lambda (),	_R ( iepdf._R() ) {}
-			void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 )
-			{
-				it_assert_debug ( A0.rows() ==mu0.length(),"" );
-				it_assert_debug ( R0.rows() ==A0.rows(),"" );
-
-				iepdf.set_parameters ( mu0,Lambda ); //
-				A = A0;
-				mu_const = mu0;
-				Re=R0;
-				Lambda = Lambda0;
-			}
-			void condition ( const vec &cond )
-			{
-				iepdf._mu() = A*cond + mu_const;
-				double zeta;
-				//ugly hack!
-				if ( ( cond.length() +1 ) ==Lambda.rows() )
-				{
-					zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) );
-				}
-				else
-				{
-					zeta = Lambda.invqform ( cond );
-				}
-				_R = Re;
-				_R*= ( 1+zeta );// / ( nu ); << nu is in Re!!!!!!
-			};
-
-	};
-	/*!
-	\brief  Gamma random walk
-
-	Mean value, \f$\mu\f$, of this density is given by \c rvc .
-	Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
-	This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
-
-	The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
-	*/
-	class mgamma : public mpdf_internal<egamma>
-	{
-		protected:
-
-			//! Constant \f$k\f$
-			double k;
-
-			//! cache of iepdf.beta
-			vec &_beta;
-
-		public:
-			//! Constructor
-			mgamma():mpdf_internal<egamma>(), k(0),
-			    _beta(iepdf._beta()) {
-			}
-
-			//! Set value of \c k
-			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", UI::compulsory);
-				UI::get(k,set,"k", UI::compulsory);
-				set_parameters(k,betatmp);
-			}
-	};
-	UIREGISTER(mgamma);
-	
-	/*!
-	\brief  Inverse-Gamma random walk
-
-	Mean value, \f$ \mu \f$, of this density is given by \c rvc .
-	Standard deviation of the random walk is proportional to one \f$ k \f$-th the mean.
-	This is achieved by setting \f$ \alpha=\mu/k^2+2 \f$ and \f$ \beta=\mu(\alpha-1)\f$.
-
-	The standard deviation of the walk is then: \f$ \mu/\sqrt(k)\f$.
-	 */
-	class migamma : public mpdf_internal<eigamma>
-	{
-		protected:
-			//! Constant \f$k\f$
-			double k;
-
-			//! cache of iepdf.alpha
-			vec &_alpha;
-
-			//! cache of iepdf.beta
-			vec &_beta;
-
-		public:
-			//! \name Constructors
-			//!@{
-			migamma():mpdf_internal<eigamma>(),
-			    k(0),
-			    _alpha(iepdf._alpha()),
-			    _beta(iepdf._beta()) {
-			}
-
-			migamma(const migamma &m):mpdf_internal<eigamma>(),
-			    k(0),
-			    _alpha(iepdf._alpha()),
-			    _beta(iepdf._beta()) {
-			}
-			//!@}
-
-			//! Set value of \c k
-			void set_parameters ( int len, double k0 )
-			{
-				k=k0;
-				iepdf.set_parameters ( ( 1.0/ ( k*k ) +2.0 ) *ones ( len ) /*alpha*/, ones ( len ) /*beta*/ );
-				dimc = dimension();
-			};
-			void condition ( const vec &val )
-			{
-				_beta=elem_mult ( val, ( _alpha-1.0 ) );
-			};
-	};
-
-
-	/*!
-	\brief  Gamma random walk around a fixed point
-
-	Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
-	\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
-
-	Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
-	This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
-
-	The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
-	*/
-	class mgamma_fix : public mgamma
-	{
-		protected:
-			//! parameter l
-			double l;
-			//! reference vector
-			vec refl;
-		public:
-			//! Constructor
-			mgamma_fix ( ) : mgamma ( ),refl () {};
-			//! Set value of \c k
-			void set_parameters ( double k0 , vec ref0, double l0 )
-			{
-				mgamma::set_parameters ( k0, ref0 );
-				refl=pow ( ref0,1.0-l0 );l=l0;
-				dimc=dimension();
-			};
-
-			void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); _beta=k/mean;};
-	};
-
-
-	/*!
-	\brief  Inverse-Gamma random walk around a fixed point
-
-	Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
-	\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]
-
-	==== Check == vv =
-	Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
-	This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
-
-	The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
-	 */
-	class migamma_ref : public migamma
-	{
-		protected:
-			//! parameter l
-			double l;
-			//! reference vector
-			vec refl;
-		public:
-			//! Constructor
-			migamma_ref ( ) : migamma (),refl ( ) {};
-			//! Set value of \c k
-			void set_parameters ( double k0 , vec ref0, double l0 )
-			{
-				migamma::set_parameters ( ref0.length(), k0 );
-				refl=pow ( ref0,1.0-l0 );
-				l=l0;
-				dimc = dimension();
-			};
-
-			void condition ( const vec &val )
-			{
-				vec mean=elem_mult ( refl,pow ( val,l ) );
-				migamma::condition ( mean );
-			};
-
-			/*! UI for migamma_ref
-
-			The migamma_ref is constructed from a structure with fields:
-			\code
-			system = {
-				type = "migamma_ref";
-				ref = [1e-5; 1e-5; 1e-2 1e-3];            // reference vector
-				l = 0.999;                                // constant l
-				k = 0.1;                                  // constant k
-				
-				// == OPTIONAL ==
-				// description of y variables
-				y = {type="rv"; names=["y", "u"];};
-				// description of u variable
-				u = {type="rv"; names=[];}
-			};
-			\endcode
-
-			Result if
-			 */
-			void from_setting( const Setting &set );
-
-			// TODO dodelat void to_setting( Setting &set ) const;
-	};
-
-
-	UIREGISTER(migamma_ref);
-
-	/*! Log-Normal probability density
-	 only allow diagonal covariances!
-
-	Density of the form \f$ \log(x)\sim \mathcal{N}(\mu,\sigma^2) \f$ , i.e.
-	\f[
-	x \sim \frac{1}{x\sigma\sqrt{2\pi}}\exp{-\frac{1}{2\sigma^2}(\log(x)-\mu)}
-	\f]
-
-	*/
-	class elognorm: public enorm<ldmat>
-	{
-		public:
-			vec sample() const {return exp ( enorm<ldmat>::sample() );};
-			vec mean() const {vec var=enorm<ldmat>::variance();return exp ( mu - 0.5*var );};
-
-	};
-
-	/*!
-	\brief  Log-Normal random walk
-
-	Mean value, \f$\mu\f$, is...
-
-	==== Check == vv =
-	Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
-	This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.
-
-	The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
-	 */
-	class mlognorm : public mpdf_internal<elognorm>
-	{
-		protected:
-			//! parameter 1/2*sigma^2
-			double sig2;
-
-			//! access
-			vec &mu;
-		public:
-			//! Constructor
-			mlognorm():mpdf_internal<elognorm>(),
-			    sig2(0),
-			    mu(iepdf._mu()) {
-			}
-
-			//! Set value of \c k
-			void set_parameters ( int size, double k )
-			{
-				sig2 = 0.5*log ( k*k+1 );
-				iepdf.set_parameters ( zeros ( size ),2*sig2*eye ( size ) );
-
-				dimc = size;
-			};
-
-			void condition ( const vec &val )
-			{
-				mu=log ( val )-sig2;//elem_mult ( refl,pow ( val,l ) );
-			};
-
-			/*! UI for mlognorm
-
-			The mlognorm is constructed from a structure with fields:
-			\code
-			system = {
-				type = "mlognorm";
-				k = 0.1;                                  // constant k
-				mu0 = [1., 1.];
-				
-				// == OPTIONAL ==
-				// description of y variables
-				y = {type="rv"; names=["y", "u"];};
-				// description of u variable
-				u = {type="rv"; names=[];}
-			};
-			\endcode
-
-			 */
-			void from_setting( const Setting &set );
-
-			// TODO dodelat void to_setting( Setting &set ) const;
-
-	};
-	
-	UIREGISTER(mlognorm);
-
-	/*! inverse Wishart density defined on Choleski decomposition
-
-	*/
-	class eWishartCh : public epdf
-	{
-		protected:
-			//! Upper-Triagle of Choleski decomposition of \f$ \Psi \f$
-			chmat Y;
-			//! dimension of matrix  \f$ \Psi \f$
-			int p;
-			//! degrees of freedom  \f$ \nu \f$
-			double delta;
-		public:
-			void set_parameters ( const mat &Y0, const double delta0 ) {Y=chmat ( Y0 );delta=delta0; p=Y.rows(); dim = p*p; }
-			mat sample_mat() const
-			{
-				mat X=zeros ( p,p );
-
-				//sample diagonal
-				for ( int i=0;i<p;i++ )
-				{
-					GamRNG.setup ( 0.5* ( delta-i ) , 0.5 ); // no +1 !! index if from 0
-#pragma omp critical
-					X ( i,i ) =sqrt ( GamRNG() );
-				}
-				//do the rest
-				for ( int i=0;i<p;i++ )
-				{
-					for ( int j=i+1;j<p;j++ )
-					{
-#pragma omp critical
-						X ( i,j ) =NorRNG.sample();
-					}
-				}
-				return X*Y._Ch();// return upper triangular part of the decomposition
-			}
-			vec sample () const
-			{
-				return vec ( sample_mat()._data(),p*p );
-			}
-			//! fast access function y0 will be copied into Y.Ch.
-			void setY ( const mat &Ch0 ) {copy_vector ( dim,Ch0._data(), Y._Ch()._data() );}
-			//! fast access function y0 will be copied into Y.Ch.
-			void _setY ( const vec &ch0 ) {copy_vector ( dim, ch0._data(), Y._Ch()._data() ); }
-			//! access function
-			const chmat& getY()const {return Y;}
-	};
-
-	class eiWishartCh: public epdf
-	{
-		protected:
-			eWishartCh W;
-			int p;
-			double delta;
-		public:
-			void set_parameters ( const mat &Y0, const double delta0) {
-				delta = delta0;
-				W.set_parameters ( inv ( Y0 ),delta0 ); 
-				dim = W.dimension(); p=Y0.rows();
-			}
-			vec sample() const {mat iCh; iCh=inv ( W.sample_mat() ); return vec ( iCh._data(),dim );}
-			void _setY ( const vec &y0 )
-			{
-				mat Ch ( p,p );
-				mat iCh ( p,p );
-				copy_vector ( dim, y0._data(), Ch._data() );
-				
-				iCh=inv ( Ch );
-				W.setY ( iCh );
-			}
-			virtual double evallog ( const vec &val ) const {
-				chmat X(p);
-				const chmat& Y=W.getY();
-				 
-				copy_vector(p*p,val._data(),X._Ch()._data());
-				chmat iX(p);X.inv(iX);
-				// compute  
-// 				\frac{ |\Psi|^{m/2}|X|^{-(m+p+1)/2}e^{-tr(\Psi X^{-1})/2} }{ 2^{mp/2}\Gamma_p(m/2)},
-				mat M=Y.to_mat()*iX.to_mat();
-				
-				double log1 = 0.5*p*(2*Y.logdet())-0.5*(delta+p+1)*(2*X.logdet())-0.5*trace(M); 
-				//Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!
-				
-/*				if (0) {
-					mat XX=X.to_mat();
-					mat YY=Y.to_mat();
-					
-					double log2 = 0.5*p*log(det(YY))-0.5*(delta+p+1)*log(det(XX))-0.5*trace(YY*inv(XX)); 
-					cout << log1 << "," << log2 << endl;
-				}*/
-				return log1;				
-			};
-			
-	};
-
-	class rwiWishartCh : public mpdf
-	{
-		protected:
-			eiWishartCh eiW;
-			//!square root of \f$ \nu-p-1 \f$ - needed for computation of \f$ \Psi \f$ from conditions
-			double sqd;
-			//reference point for diagonal
-			vec refl;
-			double l;
-			int p;
-
-		public:
-			rwiWishartCh():eiW(),
-			    sqd(0), l(0), p(0) {
-			    set_ep(eiW);
-			}
-
-			void set_parameters ( int p0, double k, vec ref0, double l0  )
-			{
-				p=p0;
-				double delta = 2/(k*k)+p+3;
-				sqd=sqrt ( delta-p-1 );
-				l=l0;
-				refl=pow(ref0,1-l);
-				
-				eiW.set_parameters ( eye ( p ),delta );
-				dimc=eiW.dimension();
-			}
-			void condition ( const vec &c ) {
-				vec z=c;
-				int ri=0;
-				for(int i=0;i<p*p;i+=(p+1)){//trace diagonal element
-					z(i) = pow(z(i),l)*refl(ri);
-					ri++;
-				}
-
-				eiW._setY ( sqd*z );
-			}
-	};
-
-//! Switch between various resampling methods.
-	enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
-	/*!
-	\brief Weighted empirical density
-
-	Used e.g. in particle filters.
-	*/
-	class eEmp: public epdf
-	{
-		protected :
-			//! Number of particles
-			int n;
-			//! Sample weights \f$w\f$
-			vec w;
-			//! Samples \f$x^{(i)}, i=1..n\f$
-			Array<vec> samples;
-		public:
-			//! \name Constructors
-			//!@{
-			eEmp ( ) :epdf ( ),w ( ),samples ( ) {};
-			//! copy constructor
-			eEmp ( const eEmp &e ) : epdf ( e ), w ( e.w ), samples ( e.samples ) {};
-			//!@}
-
-			//! Set samples and weights
-			void set_statistics ( const vec &w0, const epdf* pdf0 );
-			//! Set samples and weights
-			void set_statistics ( const epdf* pdf0 , int n ) {set_statistics ( ones ( n ) /n,pdf0 );};
-			//! Set sample
-			void set_samples ( const epdf* pdf0 );
-			//! Set sample
-			void set_parameters ( int n0, bool copy=true ) {n=n0; w.set_size ( n0,copy );samples.set_size ( n0,copy );};
-			//! Potentially dangerous, use with care.
-			vec& _w()  {return w;};
-			//! Potentially dangerous, use with care.
-			const vec& _w() const {return w;};
-			//! access function
-			Array<vec>& _samples() {return samples;};
-			//! access function
-			const Array<vec>& _samples() const {return samples;};
-			//! Function performs resampling, i.e. removal of low-weight samples and duplication of high-weight samples such that the new samples represent the same density.
-			ivec resample ( RESAMPLING_METHOD method=SYSTEMATIC );
-			//! inherited operation : NOT implemneted
-			vec sample() const {it_error ( "Not implemented" );return 0;}
-			//! inherited operation : NOT implemneted
-			double evallog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
-			vec mean() const
-			{
-				vec pom=zeros ( dim );
-				for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
-				return pom;
-			}
-			vec variance() const
-			{
-				vec pom=zeros ( dim );
-				for ( int i=0;i<n;i++ ) {pom+=pow ( samples ( i ),2 ) *w ( i );}
-				return pom-pow ( mean(),2 );
-			}
-			//! For this class, qbounds are minimum and maximum value of the population!
-			void qbounds ( vec &lb, vec &ub, double perc=0.95 ) const
-			{
-				// lb in inf so than it will be pushed below;
-				lb.set_size ( dim );
-				ub.set_size ( dim );
-				lb = std::numeric_limits<double>::infinity();
-				ub = -std::numeric_limits<double>::infinity();
-				int j;
-				for ( int i=0;i<n;i++ )
-				{
-					for ( j=0;j<dim; j++ )
-					{
-						if ( samples ( i ) ( j ) <lb ( j ) ) {lb ( j ) =samples ( i ) ( j );}
-						if ( samples ( i ) ( j ) >ub ( j ) ) {ub ( j ) =samples ( i ) ( j );}
-					}
+		//! Sample weights \f$w\f$
+		vec w;
+		//! Samples \f$x^{(i)}, i=1..n\f$
+		Array<vec> samples;
+	public:
+		//! \name Constructors
+		//!@{
+		eEmp () : epdf (), w (), samples () {};
+		//! copy constructor
+		eEmp (const eEmp &e) : epdf (e), w (e.w), samples (e.samples) {};
+		//!@}
+
+		//! Set samples and weights
+		void set_statistics (const vec &w0, const epdf &pdf0);
+		//! Set samples and weights
+		void set_statistics (const epdf &pdf0 , int n) {set_statistics (ones (n) / n, pdf0);};
+		//! Set sample
+		void set_samples (const epdf* pdf0);
+		//! Set sample
+		void set_parameters (int n0, bool copy = true) {n = n0; w.set_size (n0, copy);samples.set_size (n0, copy);};
+		//! Potentially dangerous, use with care.
+		vec& _w()  {return w;};
+		//! Potentially dangerous, use with care.
+		const vec& _w() const {return w;};
+		//! access function
+		Array<vec>& _samples() {return samples;};
+		//! access function
+		const Array<vec>& _samples() const {return samples;};
+		//! Function performs resampling, i.e. removal of low-weight samples and duplication of high-weight samples such that the new samples represent the same density.
+		ivec resample (RESAMPLING_METHOD method = SYSTEMATIC);
+		//! inherited operation : NOT implemneted
+		vec sample() const {it_error ("Not implemented");return 0;}
+		//! inherited operation : NOT implemneted
+		double evallog (const vec &val) const {it_error ("Not implemented");return 0.0;}
+		vec mean() const {
+			vec pom = zeros (dim);
+			for (int i = 0;i < n;i++) {pom += samples (i) * w (i);}
+			return pom;
+		}
+		vec variance() const {
+			vec pom = zeros (dim);
+			for (int i = 0;i < n;i++) {pom += pow (samples (i), 2) * w (i);}
+			return pom -pow (mean(), 2);
+		}
+		//! For this class, qbounds are minimum and maximum value of the population!
+		void qbounds (vec &lb, vec &ub, double perc = 0.95) const {
+			// lb in inf so than it will be pushed below;
+			lb.set_size (dim);
+			ub.set_size (dim);
+			lb = std::numeric_limits<double>::infinity();
+			ub = -std::numeric_limits<double>::infinity();
+			int j;
+			for (int i = 0;i < n;i++) {
+				for (j = 0;j < dim; j++) {
+					if (samples (i) (j) < lb (j)) {lb (j) = samples (i) (j);}
+					if (samples (i) (j) > ub (j)) {ub (j) = samples (i) (j);}
 				}
 			}
-	};
+		}
+};
 
 
 ////////////////////////
 
-	template<class sq_T>
-	void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 )
-	{
+template<class sq_T>
+void enorm<sq_T>::set_parameters (const vec &mu0, const sq_T &R0)
+{
 //Fixme test dimensions of mu0 and R0;
-		mu = mu0;
-		R = R0;
-		validate();
-	};
-
-	template<class sq_T>
-	void enorm<sq_T>::from_setting(const Setting &set){
-		epdf::from_setting(set); //reads rv
-		
-		UI::get(mu,set,"mu", UI::compulsory);
-		mat Rtmp;// necessary for conversion
-		UI::get(Rtmp,set,"R", UI::compulsory);
-		R=Rtmp; // conversion
-		validate();
-	}
-
-	template<class sq_T>
-	void enorm<sq_T>::dupdate ( mat &v, double nu )
-	{
-		//
-	};
+	mu = mu0;
+	R = R0;
+	validate();
+};
+
+template<class sq_T>
+void enorm<sq_T>::from_setting (const Setting &set)
+{
+	epdf::from_setting (set); //reads rv
+
+	UI::get (mu, set, "mu", UI::compulsory);
+	mat Rtmp;// necessary for conversion
+	UI::get (Rtmp, set, "R", UI::compulsory);
+	R = Rtmp; // conversion
+	validate();
+}
+
+template<class sq_T>
+void enorm<sq_T>::dupdate (mat &v, double nu)
+{
+	//
+};
 
 // template<class sq_T>
@@ -1204,15 +1171,15 @@
 // };
 
-	template<class sq_T>
-	vec enorm<sq_T>::sample() const
-	{
-		vec x ( dim );
+template<class sq_T>
+vec enorm<sq_T>::sample() const
+{
+	vec x (dim);
 #pragma omp critical
-		NorRNG.sample_vector ( dim,x );
-		vec smp = R.sqrt_mult ( x );
-
-		smp += mu;
-		return smp;
-	};
+	NorRNG.sample_vector (dim, x);
+	vec smp = R.sqrt_mult (x);
+
+	smp += mu;
+	return smp;
+};
 
 // template<class sq_T>
@@ -1224,19 +1191,19 @@
 // };
 
-	template<class sq_T>
-	double enorm<sq_T>::evallog_nn ( const vec &val ) const
-	{
-		// 1.83787706640935 = log(2pi)
-		double tmp=-0.5* ( R.invqform ( mu-val ) );// - lognc();
-		return  tmp;
-	};
-
-	template<class sq_T>
-	inline double enorm<sq_T>::lognc () const
-	{
-		// 1.83787706640935 = log(2pi)
-		double tmp=0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
-		return tmp;
-	};
+template<class sq_T>
+double enorm<sq_T>::evallog_nn (const vec &val) const
+{
+	// 1.83787706640935 = log(2pi)
+	double tmp = -0.5 * (R.invqform (mu - val));// - lognc();
+	return  tmp;
+};
+
+template<class sq_T>
+inline double enorm<sq_T>::lognc () const
+{
+	// 1.83787706640935 = log(2pi)
+	double tmp = 0.5 * (R.cols() * 1.83787706640935 + R.logdet());
+	return tmp;
+};
 
 
@@ -1267,67 +1234,67 @@
 
 
-	template<class sq_T>
-	enorm<sq_T>* enorm<sq_T>::marginal ( const RV &rvn ) const
-	{
-		it_assert_debug ( isnamed(), "rv description is not assigned" );
-		ivec irvn = rvn.dataind ( rv );
-
-		sq_T Rn ( R,irvn ); //select rows and columns of R
-
-		enorm<sq_T>* tmp = new enorm<sq_T>;
-		tmp->set_rv ( rvn );
-		tmp->set_parameters ( mu ( irvn ), Rn );
-		return tmp;
-	}
-
-	template<class sq_T>
-	mpdf* enorm<sq_T>::condition ( const RV &rvn ) const
-	{
-
-		it_assert_debug ( isnamed(),"rvs are not assigned" );
-
-		RV rvc = rv.subt ( rvn );
-		it_assert_debug ( ( rvc._dsize() +rvn._dsize() ==rv._dsize() ),"wrong rvn" );
-		//Permutation vector of the new R
-		ivec irvn = rvn.dataind ( rv );
-		ivec irvc = rvc.dataind ( rv );
-		ivec perm=concat ( irvn , irvc );
-		sq_T Rn ( R,perm );
-
-		//fixme - could this be done in general for all sq_T?
-		mat S=Rn.to_mat();
-		//fixme
-		int n=rvn._dsize()-1;
-		int end=R.rows()-1;
-		mat S11 = S.get ( 0,n, 0, n );
-		mat S12 = S.get ( 0, n , rvn._dsize(), end );
-		mat S22 = S.get ( rvn._dsize(), end, rvn._dsize(), end );
-
-		vec mu1 = mu ( irvn );
-		vec mu2 = mu ( irvc );
-		mat A=S12*inv ( S22 );
-		sq_T R_n ( S11 - A *S12.T() );
-
-		mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( );
-		tmp->set_rv ( rvn ); tmp->set_rvc ( rvc );
-		tmp->set_parameters ( A,mu1-A*mu2,R_n );
-		return tmp;
-	}
-
-	////
-	///////
-	template<class sq_T>				
-			void mgnorm<sq_T >::set_parameters ( fnc* g0, const sq_T &R0 ) {g=g0; this->iepdf.set_parameters ( zeros ( g->dimension() ), R0 );}
-	template<class sq_T>				
-			void mgnorm<sq_T >::condition ( const vec &cond ) {this->iepdf._mu()=g->eval ( cond );};
-
-	template<class sq_T>
-	std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_T> &ml )
-	{
-		os << "A:"<< ml.A<<endl;
-		os << "mu:"<< ml.mu_const<<endl;
-		os << "R:" << ml.iepdf->_R().to_mat() <<endl;
-		return os;
-	};
+template<class sq_T>
+enorm<sq_T>* enorm<sq_T>::marginal (const RV &rvn) const
+{
+	it_assert_debug (isnamed(), "rv description is not assigned");
+	ivec irvn = rvn.dataind (rv);
+
+	sq_T Rn (R, irvn);  //select rows and columns of R
+
+	enorm<sq_T>* tmp = new enorm<sq_T>;
+	tmp->set_rv (rvn);
+	tmp->set_parameters (mu (irvn), Rn);
+	return tmp;
+}
+
+template<class sq_T>
+mpdf* enorm<sq_T>::condition (const RV &rvn) const
+{
+
+	it_assert_debug (isnamed(), "rvs are not assigned");
+
+	RV rvc = rv.subt (rvn);
+	it_assert_debug ( (rvc._dsize() + rvn._dsize() == rv._dsize()), "wrong rvn");
+	//Permutation vector of the new R
+	ivec irvn = rvn.dataind (rv);
+	ivec irvc = rvc.dataind (rv);
+	ivec perm = concat (irvn , irvc);
+	sq_T Rn (R, perm);
+
+	//fixme - could this be done in general for all sq_T?
+	mat S = Rn.to_mat();
+	//fixme
+	int n = rvn._dsize() - 1;
+	int end = R.rows() - 1;
+	mat S11 = S.get (0, n, 0, n);
+	mat S12 = S.get (0, n , rvn._dsize(), end);
+	mat S22 = S.get (rvn._dsize(), end, rvn._dsize(), end);
+
+	vec mu1 = mu (irvn);
+	vec mu2 = mu (irvc);
+	mat A = S12 * inv (S22);
+	sq_T R_n (S11 - A *S12.T());
+
+	mlnorm<sq_T>* tmp = new mlnorm<sq_T> ();
+	tmp->set_rv (rvn); tmp->set_rvc (rvc);
+	tmp->set_parameters (A, mu1 - A*mu2, R_n);
+	return tmp;
+}
+
+////
+///////
+template<class sq_T>
+void mgnorm<sq_T >::set_parameters (fnc* g0, const sq_T &R0) {g = g0; this->iepdf.set_parameters (zeros (g->dimension()), R0);}
+template<class sq_T>
+void mgnorm<sq_T >::condition (const vec &cond) {this->iepdf._mu() = g->eval (cond);};
+
+template<class sq_T>
+std::ostream &operator<< (std::ostream &os,  mlnorm<sq_T> &ml)
+{
+	os << "A:" << ml.A << endl;
+	os << "mu:" << ml.mu_const << endl;
+	os << "R:" << ml._R() << endl;
+	return os;
+};
 
 }
Index: library/bdm/stat/merger.h
===================================================================
--- library/bdm/stat/merger.h (revision 477)
+++ library/bdm/stat/merger.h (revision 488)
@@ -178,5 +178,5 @@
 	//! Set support points from a pdf by drawing N samples
 	void set_support ( const epdf &overall, int N ) {
-		eSmp.set_statistics ( &overall, N );
+		eSmp.set_statistics ( overall, N );
 		Npoints = N;
 	}
Index: library/bdm/stat/emix.h
===================================================================
--- library/bdm/stat/emix.h (revision 487)
+++ library/bdm/stat/emix.h (revision 488)
@@ -472,38 +472,39 @@
 
 
-/*! \brief Mixture of mpdfs with constant weights, all mpdfs are of equal type
+/*! \brief Mixture of mpdfs with constant weights, all mpdfs are of equal RV and RVC
 
 */
-// class mmix : public mpdf {
-// protected:
-// 	//! Component (epdfs)
-// 	Array<mpdf*> Coms;
-// 
-// 	//!Internal epdf
-// 	emix iepdf;
-// 
-// public:
-// 	//!Default constructor
-// 	mmix() : iepdf ( ) {
-// 		set_ep ( iepdf );
-// 	}
-// 
-// 	//! Set weights \c w and components \c R
-// 	void set_parameters ( const vec &w, const Array<mpdf*> &Coms ) {
-// 		Array<epdf*> Eps ( Coms.length() );
-// 
-// 		for ( int i = 0; i < Coms.length(); i++ ) {
-// 			Eps ( i ) = Coms ( i )->e();
-// 		}
-// 
-// 		iepdf->set_parameters ( w, Eps );
-// 	}
-// 
-// 	void condition ( const vec &cond ) {
-// 		for ( int i = 0; i < Coms.length(); i++ ) {
-// 			Coms ( i )->condition ( cond );
-// 		}
-// 	};
-// };
+class mmix : public mpdf {
+protected:
+	//! Component (mpdfs)
+	Array<shared_ptr<mpdf> > Coms;
+	//!weights of the components
+	vec w;
+	//! dummy epdfs
+	epdf dummy_epdf;
+public:
+	//!Default constructor
+	mmix() : Coms(0), dummy_epdf() { set_ep(dummy_epdf);	}
+
+	//! Set weights \c w and components \c R
+	void set_parameters ( const vec &w0, const Array<shared_ptr<mpdf> > &Coms0 ) {
+		//!\TODO check if all components are OK
+		Coms = Coms0;
+		w=w0;	
+		
+		set_rv(Coms(0)->_rv());
+		set_rvc(Coms(0)->_rvc());
+	}
+	double evallogcond (const vec &dt, const vec &cond) {
+		double ll=0.0;
+		for (int i=0;i<Coms.length();i++){
+			ll+=Coms(i)->evallogcond(dt,cond);
+		}
+		return ll;
+	}
+
+	vec samplecond (const vec &cond);
+
+};
 
 }
Index: library/bdm/stat/exp_family.cpp
===================================================================
--- library/bdm/stat/exp_family.cpp (revision 487)
+++ library/bdm/stat/exp_family.cpp (revision 488)
@@ -305,15 +305,14 @@
 }
 
-void eEmp::set_statistics ( const vec &w0, const epdf* epdf0 ) {
+void eEmp::set_statistics ( const vec &w0, const epdf &epdf0 ) {
 	//it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
-	dim = epdf0->dimension();
+	dim = epdf0.dimension();
 	w = w0;
 	w /= sum ( w0 );//renormalize
 	n = w.length();
 	samples.set_size ( n );
-	dim = epdf0->dimension();
 
 	for ( int i = 0; i < n; i++ ) {
-		samples ( i ) = epdf0->sample();
+		samples ( i ) = epdf0.sample();
 	}
 }
Index: library/tests/test_particle.cpp
===================================================================
--- library/tests/test_particle.cpp (revision 477)
+++ library/tests/test_particle.cpp (revision 488)
@@ -25,5 +25,5 @@
 	euni eun;
 	eun.set_parameters ( "0", "1" );
-	emp.set_statistics ( ones ( 10 ), &eun );
+	emp.set_statistics ( ones ( 10 ), eun );
 	vec &v = emp._w();
 	Array<vec> &S = emp._samples();
Index: library/tests/test_kalman_QR.cpp
===================================================================
--- library/tests/test_kalman_QR.cpp (revision 477)
+++ library/tests/test_kalman_QR.cpp (revision 488)
@@ -69,5 +69,6 @@
 	KF_QR.set_parameters ( &evolQR, &evolQR, 100 );
 	evolQR.condition ( "1 1 1" );
-	KF_QR.set_statistics ( evolQR.e(), &KF );
+	egamma kfinit=egamma("10 10 10","1 1 1");
+	KF_QR.set_statistics ( egamma("10 10 10","1 1 1"), &KF );
 	const epdf& mpost = KF_QR.posterior();
 	const epdf& mposttr = KFtr.posterior();
Index: library/tests/testResample.cpp
===================================================================
--- library/tests/testResample.cpp (revision 477)
+++ library/tests/testResample.cpp (revision 488)
@@ -25,5 +25,5 @@
 	euni eun;
 	eun.set_parameters ( "0", "1" );
-	emp.set_statistics ( ones ( 10 ), &eun );
+	emp.set_statistics ( ones ( 10 ), eun );
 	vec &v = emp._w();
 	Array<vec> &S = emp._samples();
Index: library/tests/testSmp.cpp
===================================================================
--- library/tests/testSmp.cpp (revision 477)
+++ library/tests/testSmp.cpp (revision 488)
@@ -86,5 +86,5 @@
 	cout << "======= MMix ======== " << endl;
 	mmix mMix;
-	Array<mpdf*> mComs ( 2 );
+	Array<shared_ptr<mpdf> > mComs ( 2 );
 	mComs ( 0 ) = &mG;
 	eN.set_mu ( vec_2 ( 0.0, 0.0 ) );
@@ -94,5 +94,5 @@
 
 	Smp = mMix.samplecond_m ( mu0, N );
-	disp ( mMix.e()->mean(), zeros ( 2 ), Smp );
+	disp ( 0.5*eN.mean()+0.4*eG.mean(), zeros ( 2 ), Smp );
 
 	cout << "======= EProd ======== " << endl;
