Index: library/bdm/stat/merger.cpp
===================================================================
--- library/bdm/stat/merger.cpp (revision 713)
+++ library/bdm/stat/merger.cpp (revision 737)
@@ -5,6 +5,6 @@
 namespace bdm {
 
-merger_base::merger_base ( const Array<shared_ptr<pdf> > &S ):
-	Npoints(0), DBG(false), dbg_file(0) {
+merger_base::merger_base ( const Array<shared_ptr<pdf> > &S ) :
+		Npoints ( 0 ), DBG ( false ), dbg_file ( 0 ) {
 	set_sources ( S );
 }
@@ -103,5 +103,5 @@
 		//Re-Initialize Mixture model
 		Mix.flatten ( &Mix_init );
-		Mix.bayes_batch ( Smp_ex,empty_vec, w*Npoints );
+		Mix.bayes_batch ( Smp_ex, empty_vec, w*Npoints );
 		delete Mpred;
 		Mpred = Mix.epredictor ( ); // Allocation => must be deleted at the end!!
@@ -151,5 +151,5 @@
 			if ( pdfs ( i )->dimension() == dim ) {
 				// no need for conditioning or marginalization
-				lw_src = pdfs ( i )->evallogcond_mat ( Smp , vec(0));
+				lw_src = pdfs ( i )->evallogcond_mat ( Smp , vec ( 0 ) );
 			} else {
 				// compute likelihood of marginal on the conditional variable
@@ -187,6 +187,6 @@
 				// Compute likelihood of the partial source
 				for ( int k = 0; k < Npoints; k++ ) {
-					lw_src ( k ) += pdfs ( i )->evallogcond ( dls ( i )->pushdown ( Smp ( k ) ),  
-							 dls ( i )->get_cond ( Smp ( k ) ));
+					lw_src ( k ) += pdfs ( i )->evallogcond ( dls ( i )->pushdown ( Smp ( k ) ),
+					                dls ( i )->get_cond ( Smp ( k ) ) );
 				}
 
Index: library/bdm/stat/discrete.cpp
===================================================================
--- library/bdm/stat/discrete.cpp (revision 730)
+++ library/bdm/stat/discrete.cpp (revision 737)
@@ -122,9 +122,9 @@
 }
 
-void grid_fnc::set_values (const epdf &ep ) {
+void grid_fnc::set_values ( const epdf &ep ) {
 	if ( sup.points() > 0 ) {
-		values ( 0 ) =  exp(ep.evallog( sup.first_vec() ));
+		values ( 0 ) =  exp ( ep.evallog ( sup.first_vec() ) );
 		for ( int j = 1; j < sup.points(); j++ ) {
-			values ( j ) = exp(ep.evallog( sup.next_vec() ));
+			values ( j ) = exp ( ep.evallog ( sup.next_vec() ) );
 		}
 	}
Index: library/bdm/stat/emix.cpp
===================================================================
--- library/bdm/stat/emix.cpp (revision 693)
+++ library/bdm/stat/emix.cpp (revision 737)
@@ -38,5 +38,5 @@
 shared_ptr<epdf> emix::marginal ( const RV &rv ) const {
 	emix *tmp = new emix();
-	shared_ptr<epdf> narrow(tmp);
+	shared_ptr<epdf> narrow ( tmp );
 	marginal ( rv, *tmp );
 	return narrow;
@@ -57,5 +57,5 @@
 	bdm_assert ( isnamed(), "rvs are not assigned" );
 	mratio *tmp = new mratio ( this, rv );
-	return shared_ptr<pdf>(tmp);
+	return shared_ptr<pdf> ( tmp );
 }
 
@@ -118,5 +118,5 @@
 shared_ptr<epdf> egiwmix::marginal ( const RV &rv ) const {
 	emix *tmp = new emix();
-	shared_ptr<epdf> narrow(tmp);
+	shared_ptr<epdf> narrow ( tmp );
 	marginal ( rv, *tmp );
 	return narrow;
@@ -203,5 +203,5 @@
 };
 
-void mprod::set_elements (const Array<shared_ptr<pdf> > &mFacs ) {
+void mprod::set_elements ( const Array<shared_ptr<pdf> > &mFacs ) {
 	pdfs = mFacs;
 	dls.set_size ( mFacs.length() );
@@ -214,5 +214,5 @@
 		rvc.add ( rvx ); // add rv to common rvc
 	}
-	dimc=rvc._dsize();
+	dimc = rvc._dsize();
 
 	// rv and rvc established = > we can link them with pdfs
@@ -223,5 +223,5 @@
 }
 
-vec mmix::samplecond(const vec &cond) {
+vec mmix::samplecond ( const vec &cond ) {
 	//Sample which component
 	vec cumDist = cumsum ( w );
@@ -235,5 +235,5 @@
 	}
 
-	return Coms ( i )->samplecond(cond);
+	return Coms ( i )->samplecond ( cond );
 }
 
Index: library/bdm/stat/exp_family.h
===================================================================
--- library/bdm/stat/exp_family.h (revision 735)
+++ library/bdm/stat/exp_family.h (revision 737)
@@ -19,6 +19,5 @@
 #include "../math/chmat.h"
 
-namespace bdm
-{
+namespace bdm {
 
 
@@ -36,77 +35,79 @@
 */
 
-class eEF : public epdf
-{
-	public:
+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;
-
-		//!Evaluate normalized log-probability
-		virtual double evallog_nn (const vec &val) const {
-			bdm_error ("Not implemented");
-			return 0.0;
-		}
-
-		//!Evaluate normalized log-probability
-		virtual double evallog (const vec &val) const {
-			double tmp;
-			tmp = evallog_nn (val) - lognc();
-			return tmp;
-		}
-		//!Evaluate normalized log-probability for many samples
-		virtual vec evallog_mat (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_mat (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) {
-			bdm_error ("Not implemented");
-		}
+	//! default constructor
+	eEF () : epdf () {};
+	//! logarithm of the normalizing constant, \f$\mathcal{I}\f$
+	virtual double lognc() const = 0;
+
+	//!Evaluate normalized log-probability
+	virtual double evallog_nn ( const vec &val ) const {
+		bdm_error ( "Not implemented" );
+		return 0.0;
+	}
+
+	//!Evaluate normalized log-probability
+	virtual double evallog ( const vec &val ) const {
+		double tmp;
+		tmp = evallog_nn ( val ) - lognc();
+		return tmp;
+	}
+	//!Evaluate normalized log-probability for many samples
+	virtual vec evallog_mat ( 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_mat ( 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 ) {
+		bdm_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) {
-			bdm_error ("Not implemented");
-		}
-
-		//! Weighted update of sufficient statistics (Bayes rule)
-		virtual void bayes_weighted (const vec &data, const vec &cond=empty_vec, const double w=1.0) {};
-		//original Bayes
-		void bayes (const vec &yt, const vec &cond=empty_vec);
-
-		//!Flatten the posterior according to the given BMEF (of the same type!)
-		virtual void flatten (const BMEF * B) {
-			bdm_error ("Not implemented");
-		}
-
-		BMEF* _copy_ () const {
-			bdm_error ("function _copy_ not implemented for this BM");
-			return NULL;
-		}
+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 ) {
+		bdm_error ( "Not implemented" );
+	}
+
+	//! Weighted update of sufficient statistics (Bayes rule)
+	virtual void bayes_weighted ( const vec &data, const vec &cond = empty_vec, const double w = 1.0 ) {};
+	//original Bayes
+	void bayes ( const vec &yt, const vec &cond = empty_vec );
+
+	//!Flatten the posterior according to the given BMEF (of the same type!)
+	virtual void flatten ( const BMEF * B ) {
+		bdm_error ( "Not implemented" );
+	}
+
+	BMEF* _copy_ () const {
+		bdm_error ( "function _copy_ not implemented for this BM" );
+		return NULL;
+	}
 };
 
@@ -120,75 +121,90 @@
 */
 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);
-		/*! Create Normal density
-		\f[ f(rv) = N(\mu, R) \f]
-		from structure
-		\code
-		class = 'enorm<ldmat>', (OR) 'enorm<chmat>', (OR) 'enorm<fsqmat>';
-		mu    = [];                  // mean value
-		R     = [];                  // variance, square matrix of appropriate dimension
-		\endcode
-		*/
-		void from_setting (const Setting &root);
-		void validate() {
-			bdm_assert (mu.length() == R.rows(), "mu and R parameters do not match");
-			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());}
+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 );
+	/*! Create Normal density
+	\f[ f(rv) = N(\mu, R) \f]
+	from structure
+	\code
+	class = 'enorm<ldmat>', (OR) 'enorm<chmat>', (OR) 'enorm<fsqmat>';
+	mu    = [];                  // mean value
+	R     = [];                  // variance, square matrix of appropriate dimension
+	\endcode
+	*/
+	void from_setting ( const Setting &root );
+	void validate() {
+		bdm_assert ( mu.length() == R.rows(), "mu and R parameters do not match" );
+		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?
-		shared_ptr<pdf> condition ( const RV &rvn ) const;
-
-		// target not typed to mlnorm<sq_T, enorm<sq_T> > &
-		// because that doesn't compile (perhaps because we
-		// haven't finished defining enorm yet), but the type
-		// is required
-		void condition ( const RV &rvn, pdf &target ) const;
-
-		shared_ptr<epdf> marginal (const RV &rvn ) const;
-		void marginal ( const RV &rvn, enorm<sq_T> &target ) const;
-		//!@}
-
-		//! \name Access to attributes
-		//!@{
-
-		vec& _mu() {return mu;}
-		const vec& _mu() const {return mu;}
-		void set_mu (const vec mu0) { mu = mu0;}
-		sq_T& _R() {return R;}
-		const sq_T& _R() const {return R;}
-		//!@}
-
-};
-UIREGISTER2 (enorm, chmat);
+	shared_ptr<pdf> condition ( const RV &rvn ) const;
+
+	// target not typed to mlnorm<sq_T, enorm<sq_T> > &
+	// because that doesn't compile (perhaps because we
+	// haven't finished defining enorm yet), but the type
+	// is required
+	void condition ( const RV &rvn, pdf &target ) const;
+
+	shared_ptr<epdf> marginal ( const RV &rvn ) const;
+	void marginal ( const RV &rvn, enorm<sq_T> &target ) const;
+	//!@}
+
+	//! \name Access to attributes
+	//!@{
+
+	vec& _mu() {
+		return mu;
+	}
+	const vec& _mu() const {
+		return mu;
+	}
+	void set_mu ( const vec mu0 ) {
+		mu = mu0;
+	}
+	sq_T& _R() {
+		return R;
+	}
+	const sq_T& _R() const {
+		return R;
+	}
+	//!@}
+
+};
+UIREGISTER2 ( enorm, chmat );
 SHAREDPTR2 ( enorm, chmat );
-UIREGISTER2 (enorm, ldmat);
+UIREGISTER2 ( enorm, ldmat );
 SHAREDPTR2 ( enorm, ldmat );
-UIREGISTER2 (enorm, fsqmat);
+UIREGISTER2 ( enorm, fsqmat );
 SHAREDPTR2 ( enorm, fsqmat );
 
@@ -200,120 +216,136 @@
 *
 */
-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);
-		//!@}
-
-		vec sample() const;
-		mat sample_mat(int n) const;
-		vec mean() const;
-		vec variance() const;
-		void sample_mat(mat &Mi, chmat &Ri)const;
-			
-		void factorize(mat &M, ldmat &Vz, ldmat &Lam) const;
-		//! LS estimate of \f$\theta\f$
-		vec est_theta() const;
-
-		//! Covariance of the LS estimate
-		ldmat est_theta_cov() const;
-
-		//! expected values of the linear coefficient and the covariance matrix are written to \c M and \c R , respectively
-		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;}
-		const int & _dimx() const {return dimx;}
-			
-		/*! Create Gauss-inverse-Wishart density 
-		\f[ f(rv) = GiW(V,\nu) \f]
-		from structure
-		\code
-		class = 'egiw';
-		V     = [];               // square matrix
-		dV    = [];               // vector of diagonal of V (when V not given)
-		nu    = [];               // scalar \nu ((almost) degrees of freedom)
-								  // when missing, it will be computed to obtain proper pdf
-		dimx  = [];               // dimension of the wishart part
-		rv = RV({'name'})         // description of RV
-		rvc = RV({'name'})        // description of RV in condition
-		\endcode
-		*/
-				
-		void from_setting (const Setting &set) {
-			epdf::from_setting(set);
-			UI::get (dimx, set, "dimx", UI::compulsory);
-			if (!UI::get (nu, set, "nu", UI::optional)) {nu=-1;}
-			mat V;
-			if (!UI::get (V, set, "V", UI::optional)){
-				vec dV;
-				UI::get (dV, set, "dV", UI::compulsory);
-				set_parameters (dimx, ldmat(dV), nu);
-				
-			} else {
-				set_parameters (dimx, V, nu);
-			}
-		}
-		
-		void to_setting ( Setting& set ) const{
-			epdf::to_setting(set);
-			string s("egiw");
-			UI::save(s,set,"class");
-			UI::save(dimx,set,"dimx");
-			UI::save(V.to_mat(),set,"V");
-			UI::save(nu,set,"nu");			
-		};
-		
-		void validate(){
-			// check sizes, rvs etc.
-		}
-		void log_register( bdm::logger& L, const string& prefix ){
-			if (log_level==3){
-				root::log_register(L,prefix);
-				logrec->ids.set_length(2);
-				int th_dim=dimension()-dimx*(dimx+1)/2;
-				logrec->ids(0)=L.add_vector(RV("",th_dim), prefix + logrec->L.prefix_sep() +"mean");
-				logrec->ids(1)=L.add_vector(RV("",th_dim*th_dim),prefix + logrec->L.prefix_sep() + "variance"); 
-			} else {
-				epdf::log_register(L,prefix);
-			}
-		}
-		void log_write() const {
-			if (log_level==3){
-				mat M;
-				ldmat Lam;
-				ldmat Vz;
-				factorize(M,Vz,Lam);
-				logrec->L.log_vector(logrec->ids(0), est_theta() );
-				logrec->L.log_vector(logrec->ids(1), cvectorize(est_theta_cov().to_mat()));
-			} else {
-				epdf::log_write();
-			}
-			
-		}
-		//!@}
+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 );
+	//!@}
+
+	vec sample() const;
+	mat sample_mat ( int n ) const;
+	vec mean() const;
+	vec variance() const;
+	void sample_mat ( mat &Mi, chmat &Ri ) const;
+
+	void factorize ( mat &M, ldmat &Vz, ldmat &Lam ) const;
+	//! LS estimate of \f$\theta\f$
+	vec est_theta() const;
+
+	//! Covariance of the LS estimate
+	ldmat est_theta_cov() const;
+
+	//! expected values of the linear coefficient and the covariance matrix are written to \c M and \c R , respectively
+	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;
+	}
+	const int & _dimx() const {
+		return dimx;
+	}
+
+	/*! Create Gauss-inverse-Wishart density
+	\f[ f(rv) = GiW(V,\nu) \f]
+	from structure
+	\code
+	class = 'egiw';
+	V     = [];               // square matrix
+	dV    = [];               // vector of diagonal of V (when V not given)
+	nu    = [];               // scalar \nu ((almost) degrees of freedom)
+							  // when missing, it will be computed to obtain proper pdf
+	dimx  = [];               // dimension of the wishart part
+	rv = RV({'name'})         // description of RV
+	rvc = RV({'name'})        // description of RV in condition
+	\endcode
+	*/
+
+	void from_setting ( const Setting &set ) {
+		epdf::from_setting ( set );
+		UI::get ( dimx, set, "dimx", UI::compulsory );
+		if ( !UI::get ( nu, set, "nu", UI::optional ) ) {
+			nu = -1;
+		}
+		mat V;
+		if ( !UI::get ( V, set, "V", UI::optional ) ) {
+			vec dV;
+			UI::get ( dV, set, "dV", UI::compulsory );
+			set_parameters ( dimx, ldmat ( dV ), nu );
+
+		} else {
+			set_parameters ( dimx, V, nu );
+		}
+	}
+
+	void to_setting ( Setting& set ) const {
+		epdf::to_setting ( set );
+		string s ( "egiw" );
+		UI::save ( s, set, "class" );
+		UI::save ( dimx, set, "dimx" );
+		UI::save ( V.to_mat(), set, "V" );
+		UI::save ( nu, set, "nu" );
+	};
+
+	void validate() {
+		// check sizes, rvs etc.
+	}
+	void log_register ( bdm::logger& L, const string& prefix ) {
+		if ( log_level == 3 ) {
+			root::log_register ( L, prefix );
+			logrec->ids.set_length ( 2 );
+			int th_dim = dimension() - dimx * ( dimx + 1 ) / 2;
+			logrec->ids ( 0 ) = L.add_vector ( RV ( "", th_dim ), prefix + logrec->L.prefix_sep() + "mean" );
+			logrec->ids ( 1 ) = L.add_vector ( RV ( "", th_dim * th_dim ), prefix + logrec->L.prefix_sep() + "variance" );
+		} else {
+			epdf::log_register ( L, prefix );
+		}
+	}
+	void log_write() const {
+		if ( log_level == 3 ) {
+			mat M;
+			ldmat Lam;
+			ldmat Vz;
+			factorize ( M, Vz, Lam );
+			logrec->L.log_vector ( logrec->ids ( 0 ), est_theta() );
+			logrec->L.log_vector ( logrec->ids ( 1 ), cvectorize ( est_theta_cov().to_mat() ) );
+		} else {
+			epdf::log_write();
+		}
+
+	}
+	//!@}
 };
 UIREGISTER ( egiw );
@@ -328,77 +360,90 @@
 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();
-		}
-		//!@}
-
-		//! using sampling procedure from wikipedia
-		vec sample() const {
-			vec y(beta.length());
-			for (int i=0; i<beta.length(); i++){
-				GamRNG.setup(beta(i),1);
-				#pragma omp critical
-				y(i)=GamRNG();
-			}
-			return y/sum(y);
-		}
-
-		vec mean() const {return beta / sum (beta);};
-		vec variance() const {double gamma = sum (beta); return elem_mult (beta, (gamma-beta)) / (gamma*gamma* (gamma + 1));}
-		//! In this instance, val is ...
-		double evallog_nn (const vec &val) const {
-			double tmp; tmp = (beta - 1) * log (val);
-			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);
-			return tmp;
-		}
-
-		//!access function
-		vec& _beta()  {return beta;}
-		/*! configuration structure 
-		\code 
-		class = 'eDirich';    
-		beta  = [];           //parametr beta
-		\endcode
-		*/
-		void from_setting(const Setting &set){
-			epdf::from_setting(set);
-			UI::get(beta,set, "beta", UI::compulsory);
-			validate();
-		}
-		void validate() {
-			//check rv
-			dim = beta.length();
-		}
-};
-UIREGISTER(eDirich);
+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();
+	}
+	//!@}
+
+	//! using sampling procedure from wikipedia
+	vec sample() const {
+		vec y ( beta.length() );
+		for ( int i = 0; i < beta.length(); i++ ) {
+			GamRNG.setup ( beta ( i ), 1 );
+#pragma omp critical
+			y ( i ) = GamRNG();
+		}
+		return y / sum ( y );
+	}
+
+	vec mean() const {
+		return beta / sum ( beta );
+	};
+	vec variance() const {
+		double gamma = sum ( beta );
+		return elem_mult ( beta, ( gamma - beta ) ) / ( gamma*gamma* ( gamma + 1 ) );
+	}
+	//! In this instance, val is ...
+	double evallog_nn ( const vec &val ) const {
+		double tmp;
+		tmp = ( beta - 1 ) * log ( val );
+		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 );
+		return tmp;
+	}
+
+	//!access function
+	vec& _beta()  {
+		return beta;
+	}
+	/*! configuration structure
+	\code
+	class = 'eDirich';
+	beta  = [];           //parametr beta
+	\endcode
+	*/
+	void from_setting ( const Setting &set ) {
+		epdf::from_setting ( set );
+		UI::get ( beta, set, "beta", UI::compulsory );
+		validate();
+	}
+	void validate() {
+		//check rv
+		dim = beta.length();
+	}
+};
+UIREGISTER ( eDirich );
 
 /*! Random Walk on Dirichlet
-Using simple assignment 
+Using simple assignment
  \f[ \beta = rvc / k + \beta_c \f]
  hence, mean value = rvc, variance = (k+1)*mean*mean;
- 
+
  The greater k is, the greater is the variance of the random walk;
- 
+
  \f$ \beta_c \f$ is used as regularizing element to avoid corner cases, i.e. when one element of rvc is zero.
  By default is it set to 0.1;
@@ -406,110 +451,131 @@
 
 class mDirich: public pdf_internal<eDirich> {
-	protected:
-		//! constant \f$ k \f$ of the random walk
-		double k;
-		//! cache of beta_i
-		vec &_beta;
-		//! stabilizing coefficient \f$ \beta_c \f$
-		vec betac;
-	public:
-		mDirich(): pdf_internal<eDirich>(), _beta(iepdf._beta()){};
-		void condition (const vec &val) {_beta =  val/k+betac; };
-		/*! Create Dirichlet random walk 
-		\f[ f(rv|rvc) = Di(rvc*k) \f]
-		from structure
-		\code
-		class = 'mDirich';
-		k = 1;                      // multiplicative constant k
-		--- optional ---
-		rv = RV({'name'},size)      // description of RV
-		beta0 = [];                 // initial value of beta
-		betac = [];                 // initial value of beta
-		\endcode
-		*/
-		void from_setting (const Setting &set) {
-			pdf::from_setting (set); // reads rv and rvc
-			if (_rv()._dsize()>0){
-				rvc = _rv().copy_t(-1);
-			}
-			vec beta0; 
-			if (!UI::get (beta0, set, "beta0", UI::optional)){
-				beta0 = ones(_rv()._dsize());
-			}
-			if (!UI::get (betac, set, "betac", UI::optional)){
-				betac = 0.1*ones(_rv()._dsize());
-			}
-			_beta = beta0;
-			
-			UI::get (k, set, "k", UI::compulsory);
-			validate();
-		}
-		void validate() { 
-			pdf_internal<eDirich>::validate();
-			bdm_assert(_beta.length()==betac.length(),"beta0 and betac are not compatible");
-			if (_rv()._dsize()>0){
-				bdm_assert( (_rv()._dsize()==dimension()) , "Size of rv does not match with beta");
-			}
-			dimc = _beta.length();
-		};
-};
-UIREGISTER(mDirich);
+protected:
+	//! constant \f$ k \f$ of the random walk
+	double k;
+	//! cache of beta_i
+	vec &_beta;
+	//! stabilizing coefficient \f$ \beta_c \f$
+	vec betac;
+public:
+	mDirich() : pdf_internal<eDirich>(), _beta ( iepdf._beta() ) {};
+	void condition ( const vec &val ) {
+		_beta =  val / k + betac;
+	};
+	/*! Create Dirichlet random walk
+	\f[ f(rv|rvc) = Di(rvc*k) \f]
+	from structure
+	\code
+	class = 'mDirich';
+	k = 1;                      // multiplicative constant k
+	--- optional ---
+	rv = RV({'name'},size)      // description of RV
+	beta0 = [];                 // initial value of beta
+	betac = [];                 // initial value of beta
+	\endcode
+	*/
+	void from_setting ( const Setting &set ) {
+		pdf::from_setting ( set ); // reads rv and rvc
+		if ( _rv()._dsize() > 0 ) {
+			rvc = _rv().copy_t ( -1 );
+		}
+		vec beta0;
+		if ( !UI::get ( beta0, set, "beta0", UI::optional ) ) {
+			beta0 = ones ( _rv()._dsize() );
+		}
+		if ( !UI::get ( betac, set, "betac", UI::optional ) ) {
+			betac = 0.1 * ones ( _rv()._dsize() );
+		}
+		_beta = beta0;
+
+		UI::get ( k, set, "k", UI::compulsory );
+		validate();
+	}
+	void validate() {
+		pdf_internal<eDirich>::validate();
+		bdm_assert ( _beta.length() == betac.length(), "beta0 and betac are not compatible" );
+		if ( _rv()._dsize() > 0 ) {
+			bdm_assert ( ( _rv()._dsize() == dimension() ) , "Size of rv does not match with beta" );
+		}
+		dimc = _beta.length();
+	};
+};
+UIREGISTER ( mDirich );
 
 //! \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 &yt, const vec &cond=empty_vec) {
-			if (frg < 1.0) {beta *= frg;last_lognc = est.lognc();}
-			beta += yt;
-			if (evalll) {ll = est.lognc() - last_lognc;}
-		}
-		double logpred (const vec &yt) 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 += yt;
-			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();}
-		}
-		//! return correctly typed posterior (covariant return)
-		const eDirich& posterior() const {return est;};
-		//! constructor function		
-		void set_parameters (const vec &beta0) {
-			est.set_parameters (beta0);
-			if (evalll) {last_lognc = est.lognc();}
-		}
-		void to_setting(Setting &set) const{
-			BMEF::to_setting(set);
-			Setting& prior= set.add("prior", Setting::TypeGroup);
-			est.to_setting(prior);
-		}
+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 &yt, const vec &cond = empty_vec ) {
+		if ( frg < 1.0 ) {
+			beta *= frg;
+			last_lognc = est.lognc();
+		}
+		beta += yt;
+		if ( evalll ) {
+			ll = est.lognc() - last_lognc;
+		}
+	}
+	double logpred ( const vec &yt ) 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 += yt;
+		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();
+		}
+	}
+	//! return correctly typed posterior (covariant return)
+	const eDirich& posterior() const {
+		return est;
+	};
+	//! constructor function
+	void set_parameters ( const vec &beta0 ) {
+		est.set_parameters ( beta0 );
+		if ( evalll ) {
+			last_lognc = est.lognc();
+		}
+	}
+	void to_setting ( Setting &set ) const {
+		BMEF::to_setting ( set );
+		Setting& prior = set.add ( "prior", Setting::TypeGroup );
+		est.to_setting ( prior );
+	}
 };
 
@@ -523,51 +589,63 @@
 */
 
-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;
-		double evallog (const vec &val) const;
-		double lognc () const;
-		//! Returns pointer to internal alpha. Potentially dengerous: use with care!
-		vec& _alpha() {return alpha;}
-		//! Returns pointer to internal beta. Potentially dengerous: use with care!
-		vec& _beta() {return beta;}
-		vec mean() const {return elem_div (alpha, beta);}
-		vec variance() const {return elem_div (alpha, elem_mult (beta, beta)); }
-
-		/*! Create Gamma density 
-		\f[ f(rv|rvc) = \Gamma(\alpha, \beta) \f]
-		from structure
-		\code
-		class = 'egamma';
-		alpha = [...];         // vector of alpha
-		beta = [...];          // vector of beta
-		rv = RV({'name'})      // description of RV
-		\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() {
-			bdm_assert (alpha.length() == beta.length(), "parameters do not match");
-			dim = alpha.length();
-		}
-};
-UIREGISTER (egamma);
+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;
+	double evallog ( const vec &val ) const;
+	double lognc () const;
+	//! Returns pointer to internal alpha. Potentially dengerous: use with care!
+	vec& _alpha() {
+		return alpha;
+	}
+	//! Returns pointer to internal beta. Potentially dengerous: use with care!
+	vec& _beta() {
+		return beta;
+	}
+	vec mean() const {
+		return elem_div ( alpha, beta );
+	}
+	vec variance() const {
+		return elem_div ( alpha, elem_mult ( beta, beta ) );
+	}
+
+	/*! Create Gamma density
+	\f[ f(rv|rvc) = \Gamma(\alpha, \beta) \f]
+	from structure
+	\code
+	class = 'egamma';
+	alpha = [...];         // vector of alpha
+	beta = [...];          // vector of beta
+	rv = RV({'name'})      // description of RV
+	\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() {
+		bdm_assert ( alpha.length() == beta.length(), "parameters do not match" );
+		dim = alpha.length();
+	}
+};
+UIREGISTER ( egamma );
 SHAREDPTR ( egamma );
 
@@ -588,17 +666,23 @@
  */
 
-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);}
+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 );
+	}
 };
 /*
@@ -619,94 +703,100 @@
 //!  Uniform distributed density on a rectangular support
 
-class euni: public epdf
-{
-	protected:
+class euni: public epdf {
+protected:
 //! lower bound on support
-		vec low;
+	vec low;
 //! upper bound on support
-		vec high;
+	vec high;
 //! internal
-		vec distance;
+	vec distance;
 //! normalizing coefficients
-		double nk;
+	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;
-			low = low0;
-			high = high0;
-			nk = prod (1.0 / distance);
-			lnk = log (nk);
-			dim = low.length();
-		}
-		//!@}
-
-		double evallog (const vec &val) const  {
-			if (any (val < low) && any (val > high)) {return -inf;}
-			else return lnk;
-		}
-		vec sample() const {
-			vec smp (dim);
+	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;
+		low = low0;
+		high = high0;
+		nk = prod ( 1.0 / distance );
+		lnk = log ( nk );
+		dim = low.length();
+	}
+	//!@}
+
+	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;}
-		/*! Create Uniform density 
-		\f[ f(rv) = U(low,high) \f]
-		from structure
-		 \code
-		 class = 'euni'
-		 high = [...];          // vector of upper bounds
-		 low = [...];           // vector of lower bounds
-		 rv = RV({'name'});     // 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);
-			set_parameters(low,high);
-			validate();
-		}
-		void validate() {
-			bdm_assert(high.length()==low.length(), "Incompatible high and low vectors");
-			dim = high.length();
-			bdm_assert (min (distance) > 0.0, "bad support");
-		}
-};
-UIREGISTER(euni);
+		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;
+	}
+	/*! Create Uniform density
+	\f[ f(rv) = U(low,high) \f]
+	from structure
+	 \code
+	 class = 'euni'
+	 high = [...];          // vector of upper bounds
+	 low = [...];           // vector of lower bounds
+	 rv = RV({'name'});     // 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 );
+		set_parameters ( low, high );
+		validate();
+	}
+	void validate() {
+		bdm_assert ( high.length() == low.length(), "Incompatible high and low vectors" );
+		dim = high.length();
+		bdm_assert ( min ( distance ) > 0.0, "bad support" );
+	}
+};
+UIREGISTER ( euni );
 
 //! Uniform density with conditional mean value
-class mguni : public pdf_internal<euni>{
+class mguni : public pdf_internal<euni> {
 	//! function of the mean value
 	shared_ptr<fnc> mean;
 	//! distance from mean to both sides
 	vec delta;
-	public:
-		void condition(const vec &cond){
-			vec mea=mean->eval(cond);
-			iepdf.set_parameters(mea-delta,mea+delta);
-		}
-		//! load from
-		void from_setting(const Setting &set){
-			pdf::from_setting(set); //reads rv and rvc
-			UI::get(delta,set,"delta",UI::compulsory);
-			mean = UI::build<fnc>(set,"mean",UI::compulsory);
-			
-			iepdf.set_parameters(-delta,delta);
-			dimc = mean->dimensionc();
-			validate();
-		}
-};
-UIREGISTER(mguni);
+public:
+	void condition ( const vec &cond ) {
+		vec mea = mean->eval ( cond );
+		iepdf.set_parameters ( mea - delta, mea + delta );
+	}
+	//! load from
+	void from_setting ( const Setting &set ) {
+		pdf::from_setting ( set ); //reads rv and rvc
+		UI::get ( delta, set, "delta", UI::compulsory );
+		mean = UI::build<fnc> ( set, "mean", UI::compulsory );
+
+		iepdf.set_parameters ( -delta, delta );
+		dimc = mean->dimensionc();
+		validate();
+	}
+};
+UIREGISTER ( mguni );
 /*!
  \brief Normal distributed linear function with linear function of mean value;
@@ -715,135 +805,141 @@
 */
 template < class sq_T, template <typename> class TEpdf = enorm >
-class mlnorm : public pdf_internal< TEpdf<sq_T> >
-{
-	protected:
-		//! Internal epdf that arise by conditioning on \c rvc
-		mat A;
-		//! Constant additive term 
-		vec mu_const;
+class mlnorm : public pdf_internal< TEpdf<sq_T> > {
+protected:
+	//! Internal epdf that arise by conditioning on \c rvc
+	mat A;
+	//! Constant additive term
+	vec mu_const;
 //			vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
-	public:
-		//! \name Constructors
-		//!@{
-		mlnorm() : pdf_internal< TEpdf<sq_T> >() {};
-		mlnorm (const mat &A, const vec &mu0, const sq_T &R) : pdf_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) {	
-			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;
+public:
+	//! \name Constructors
+	//!@{
+	mlnorm() : pdf_internal< TEpdf<sq_T> >() {};
+	mlnorm ( const mat &A, const vec &mu0, const sq_T &R ) : pdf_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 ) {
+		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
-		const vec& _mu_const() const {return mu_const;}
-		//!access function
-		const mat& _A() const {return A;}
-		//!access function
-		mat _R() const { return this->iepdf._R().to_mat(); }
-		//!access function
-		sq_T __R() const { return this->iepdf._R(); }
-		
-		//! Debug stream
-		template<typename sq_M>
-		friend std::ostream &operator<< (std::ostream &os,  mlnorm<sq_M, enorm> &ml);
-
-		/*! Create Normal density with linear function of mean value 
-		 \f[ f(rv|rvc) = N(A*rvc+const, R) \f]
-		from structure
-		\code
-		class = 'mlnorm<ldmat>', (OR) 'mlnorm<chmat>', (OR) 'mlnorm<fsqmat>';
-		A     = [];                  // matrix or vector of appropriate dimension
-		const = [];                  // vector of constant term
-		R     = [];                  // square matrix of appropriate dimension
-		\endcode
-		*/
-		void from_setting (const Setting &set) {
-			pdf::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);
-			validate();
-		};
-		void validate() {
-			pdf_internal<TEpdf<sq_T> >::validate();
-			bdm_assert (A.rows() == mu_const.length(), "mlnorm: A vs. mu mismatch");
-			bdm_assert (A.rows() == _R().rows(), "mlnorm: A vs. R mismatch");
-			
-		}
-};
-UIREGISTER2 (mlnorm,ldmat);
+	}
+
+	//!access function
+	const vec& _mu_const() const {
+		return mu_const;
+	}
+	//!access function
+	const mat& _A() const {
+		return A;
+	}
+	//!access function
+	mat _R() const {
+		return this->iepdf._R().to_mat();
+	}
+	//!access function
+	sq_T __R() const {
+		return this->iepdf._R();
+	}
+
+	//! Debug stream
+	template<typename sq_M>
+	friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M, enorm> &ml );
+
+	/*! Create Normal density with linear function of mean value
+	 \f[ f(rv|rvc) = N(A*rvc+const, R) \f]
+	from structure
+	\code
+	class = 'mlnorm<ldmat>', (OR) 'mlnorm<chmat>', (OR) 'mlnorm<fsqmat>';
+	A     = [];                  // matrix or vector of appropriate dimension
+	const = [];                  // vector of constant term
+	R     = [];                  // square matrix of appropriate dimension
+	\endcode
+	*/
+	void from_setting ( const Setting &set ) {
+		pdf::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 );
+		validate();
+	};
+	void validate() {
+		pdf_internal<TEpdf<sq_T> >::validate();
+		bdm_assert ( A.rows() == mu_const.length(), "mlnorm: A vs. mu mismatch" );
+		bdm_assert ( A.rows() == _R().rows(), "mlnorm: A vs. R mismatch" );
+
+	}
+};
+UIREGISTER2 ( mlnorm, ldmat );
 SHAREDPTR2 ( mlnorm, ldmat );
-UIREGISTER2 (mlnorm,fsqmat);
+UIREGISTER2 ( mlnorm, fsqmat );
 SHAREDPTR2 ( mlnorm, fsqmat );
-UIREGISTER2 (mlnorm, chmat);
+UIREGISTER2 ( mlnorm, chmat );
 SHAREDPTR2 ( mlnorm, chmat );
 
 //! pdf with general function for mean value
 template<class sq_T>
-class mgnorm : public pdf_internal< enorm< sq_T > >
-{
-	private:
+class mgnorm : public pdf_internal< enorm< sq_T > > {
+private:
 //			vec &mu; WHY NOT?
-		shared_ptr<fnc> g;
-
-	public:
-		//!default constructor
-		mgnorm() : pdf_internal<enorm<sq_T> >() { }
-		//!set mean function
-		inline void set_parameters (const shared_ptr<fnc> &g0, const sq_T &R0);
-		inline void condition (const vec &cond);
-
-
-		/*! Create Normal density with given function of mean value 
-		\f[ f(rv|rvc) = N( g(rvc), R) \f]
-		from structure
-		\code
-		class = 'mgnorm';
-		g.class =  'fnc';      // function for mean value evolution
-		g._fields_of_fnc = ...;
-
-		R = [1, 0;             // covariance matrix
-			0, 1];
-			--OR --
-		dR = [1, 1];           // diagonal of cavariance matrix
-
-		rv = RV({'name'})      // description of RV
-		rvc = RV({'name'})     // description of RV in condition
-		\endcode
-		*/
-
-		void from_setting (const Setting &set) {
-			pdf::from_setting(set);
-			shared_ptr<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);
-			validate();
-		}
-		void validate() {
-			bdm_assert(g->dimension()==this->dimension(),"incompatible function");
-		}
-};
-
-UIREGISTER2 (mgnorm, chmat);
+	shared_ptr<fnc> g;
+
+public:
+	//!default constructor
+	mgnorm() : pdf_internal<enorm<sq_T> >() { }
+	//!set mean function
+	inline void set_parameters ( const shared_ptr<fnc> &g0, const sq_T &R0 );
+	inline void condition ( const vec &cond );
+
+
+	/*! Create Normal density with given function of mean value
+	\f[ f(rv|rvc) = N( g(rvc), R) \f]
+	from structure
+	\code
+	class = 'mgnorm';
+	g.class =  'fnc';      // function for mean value evolution
+	g._fields_of_fnc = ...;
+
+	R = [1, 0;             // covariance matrix
+		0, 1];
+		--OR --
+	dR = [1, 1];           // diagonal of cavariance matrix
+
+	rv = RV({'name'})      // description of RV
+	rvc = RV({'name'})     // description of RV in condition
+	\endcode
+	*/
+
+	void from_setting ( const Setting &set ) {
+		pdf::from_setting ( set );
+		shared_ptr<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 );
+		validate();
+	}
+	void validate() {
+		bdm_assert ( g->dimension() == this->dimension(), "incompatible function" );
+	}
+};
+
+UIREGISTER2 ( mgnorm, chmat );
 SHAREDPTR2 ( mgnorm, chmat );
 
@@ -856,46 +952,45 @@
 Perhaps a moment-matching technique?
 */
-class mlstudent : public mlnorm<ldmat, enorm>
-{
-	protected:
-		//! Variable \f$ \Lambda \f$ from theory
-		ldmat Lambda;
-		//! Reference to variable \f$ R \f$
-		ldmat &_R;
-		//! Variable \f$ R_e \f$
-		ldmat Re;
-	public:
-		mlstudent () : mlnorm<ldmat, enorm> (),
-				Lambda (),	_R (iepdf._R()) {}
-		//! constructor function
-		void set_parameters (const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0) {
-			iepdf.set_parameters (mu0, R0);// was Lambda, why?
-			A = A0;
-			mu_const = mu0;
-			Re = R0;
-			Lambda = Lambda0;
-		}
-		void condition (const vec &cond) {
-			if (cond.length()>0) {
-				iepdf._mu() = A * cond + mu_const;
-			} else {
-				iepdf._mu() =  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!!!!!!
-		};
-
-		void validate() {
-			bdm_assert (A.rows() == mu_const.length(), "mlstudent: A vs. mu mismatch");
-			bdm_assert (_R.rows() == A.rows(), "mlstudent: A vs. R mismatch");
-			
-		}
+class mlstudent : public mlnorm<ldmat, enorm> {
+protected:
+	//! Variable \f$ \Lambda \f$ from theory
+	ldmat Lambda;
+	//! Reference to variable \f$ R \f$
+	ldmat &_R;
+	//! Variable \f$ R_e \f$
+	ldmat Re;
+public:
+	mlstudent () : mlnorm<ldmat, enorm> (),
+			Lambda (),	_R ( iepdf._R() ) {}
+	//! constructor function
+	void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 ) {
+		iepdf.set_parameters ( mu0, R0 );// was Lambda, why?
+		A = A0;
+		mu_const = mu0;
+		Re = R0;
+		Lambda = Lambda0;
+	}
+	void condition ( const vec &cond ) {
+		if ( cond.length() > 0 ) {
+			iepdf._mu() = A * cond + mu_const;
+		} else {
+			iepdf._mu() =  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!!!!!!
+	};
+
+	void validate() {
+		bdm_assert ( A.rows() == mu_const.length(), "mlstudent: A vs. mu mismatch" );
+		bdm_assert ( _R.rows() == A.rows(), "mlstudent: A vs. R mismatch" );
+
+	}
 };
 /*!
@@ -908,52 +1003,53 @@
 The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
 */
-class mgamma : public pdf_internal<egamma>
-{
-	protected:
-
-		//! Constant \f$k\f$
-		double k;
-
-		//! cache of iepdf.beta
-		vec &_beta;
-
-	public:
-		//! Constructor
-		mgamma() : pdf_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;};
-		/*! Create Gamma density with conditional mean value 
-		\f[ f(rv|rvc) = \Gamma(k, k/rvc) \f]
-		from structure
-		\code
-		  class = 'mgamma';
-		  beta = [...];          // vector of initial alpha
-		  k = 1.1;               // multiplicative constant k
-		  rv = RV({'name'})      // description of RV
-		  rvc = RV({'name'})     // description of RV in condition
-		 \endcode
-		*/
-		void from_setting (const Setting &set) {
-			pdf::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);
-			validate();
-		}
-		void validate() {
-			pdf_internal<egamma>::validate();
-			
-			dim = _beta.length();
-			dimc = _beta.length();
-		}
-};
-UIREGISTER (mgamma);
-SHAREDPTR (mgamma);
+class mgamma : public pdf_internal<egamma> {
+protected:
+
+	//! Constant \f$k\f$
+	double k;
+
+	//! cache of iepdf.beta
+	vec &_beta;
+
+public:
+	//! Constructor
+	mgamma() : pdf_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;
+	};
+	/*! Create Gamma density with conditional mean value
+	\f[ f(rv|rvc) = \Gamma(k, k/rvc) \f]
+	from structure
+	\code
+	  class = 'mgamma';
+	  beta = [...];          // vector of initial alpha
+	  k = 1.1;               // multiplicative constant k
+	  rv = RV({'name'})      // description of RV
+	  rvc = RV({'name'})     // description of RV in condition
+	 \endcode
+	*/
+	void from_setting ( const Setting &set ) {
+		pdf::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 );
+		validate();
+	}
+	void validate() {
+		pdf_internal<egamma>::validate();
+
+		dim = _beta.length();
+		dimc = _beta.length();
+	}
+};
+UIREGISTER ( mgamma );
+SHAREDPTR ( mgamma );
 
 /*!
@@ -966,41 +1062,40 @@
 The standard deviation of the walk is then: \f$ \mu/\sqrt(k)\f$.
  */
-class migamma : public pdf_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() : pdf_internal<eigamma>(),
-				k (0),
-				_alpha (iepdf._alpha()),
-				_beta (iepdf._beta()) {
-		}
-
-		migamma (const migamma &m) : pdf_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));
-		};
+class migamma : public pdf_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() : pdf_internal<eigamma>(),
+			k ( 0 ),
+			_alpha ( iepdf._alpha() ),
+			_beta ( iepdf._beta() ) {
+	}
+
+	migamma ( const migamma &m ) : pdf_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 ) );
+	};
 };
 
@@ -1017,22 +1112,25 @@
 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;};
+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;
+	};
 };
 
@@ -1050,48 +1148,47 @@
 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);
-		};
-
-
-		/*! Create inverse-Gamma density with conditional mean value 
-		\f[ f(rv|rvc) = i\Gamma(k, k/(rvc^l \circ ref^{(1-l)}) \f]
-		from structure
-		\code
-		class = 'migamma_ref';
-		ref = [1e-5; 1e-5; 1e-2 1e-3];            // reference vector
-		l = 0.999;                                // constant l
-		k = 0.1;                                  // constant k
-		rv = RV({'name'})                         // description of RV
-		rvc = RV({'name'})                        // description of RV in condition
-		\endcode
-		 */
-		void from_setting (const Setting &set);
-
-		// TODO dodelat void to_setting( Setting &set ) const;
-};
-
-
-UIREGISTER (migamma_ref);
-SHAREDPTR (migamma_ref);
+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 );
+	};
+
+
+	/*! Create inverse-Gamma density with conditional mean value
+	\f[ f(rv|rvc) = i\Gamma(k, k/(rvc^l \circ ref^{(1-l)}) \f]
+	from structure
+	\code
+	class = 'migamma_ref';
+	ref = [1e-5; 1e-5; 1e-2 1e-3];            // reference vector
+	l = 0.999;                                // constant l
+	k = 0.1;                                  // constant k
+	rv = RV({'name'})                         // description of RV
+	rvc = RV({'name'})                        // description of RV in condition
+	\endcode
+	 */
+	void from_setting ( const Setting &set );
+
+	// TODO dodelat void to_setting( Setting &set ) const;
+};
+
+
+UIREGISTER ( migamma_ref );
+SHAREDPTR ( migamma_ref );
 
 /*! Log-Normal probability density
@@ -1105,9 +1202,13 @@
 Function from_setting loads mu and R in the same way as it does for enorm<>!
 */
-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);};
+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 );
+	};
 
 };
@@ -1119,96 +1220,110 @@
 
  */
-class mlognorm : public pdf_internal<elognorm>
-{
-	protected:
-		//! parameter 1/2*sigma^2
-		double sig2;
-
-		//! access
-		vec &mu;
-	public:
-		//! Constructor
-		mlognorm() : pdf_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 ) );
-		};
-
-		/*! Create logNormal random Walk
-		\f[ f(rv|rvc) = log\mathcal{N}( \log(rvc)-0.5\log(k^2+1), k I) \f]
-		from structure
-		\code
-		class = 'mlognorm';
-		k   = 0.1;               // "variance" k
-		mu0 = 0.1;               // Initial value of mean
-		rv  = RV({'name'})       // description of RV
-		rvc = RV({'name'})       // description of RV in condition
-		\endcode
-		*/
-		void from_setting (const Setting &set);
-
-		// TODO dodelat void to_setting( Setting &set ) const;
-
-};
-
-UIREGISTER (mlognorm);
-SHAREDPTR (mlognorm);
+class mlognorm : public pdf_internal<elognorm> {
+protected:
+	//! parameter 1/2*sigma^2
+	double sig2;
+
+	//! access
+	vec &mu;
+public:
+	//! Constructor
+	mlognorm() : pdf_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 ) );
+	};
+
+	/*! Create logNormal random Walk
+	\f[ f(rv|rvc) = log\mathcal{N}( \log(rvc)-0.5\log(k^2+1), k I) \f]
+	from structure
+	\code
+	class = 'mlognorm';
+	k   = 0.1;               // "variance" k
+	mu0 = 0.1;               // Initial value of mean
+	rv  = RV({'name'})       // description of RV
+	rvc = RV({'name'})       // description of RV in condition
+	\endcode
+	*/
+	void from_setting ( const Setting &set );
+
+	// TODO dodelat void to_setting( Setting &set ) const;
+
+};
+
+UIREGISTER ( mlognorm );
+SHAREDPTR ( 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:
-		//! Set internal structures
-		void set_parameters (const mat &Y0, const double delta0) {Y = chmat (Y0);delta = delta0; p = Y.rows(); dim = p * p; }
-		//! Set internal structures
-		void set_parameters (const chmat &Y0, const double delta0) {Y = Y0;delta = delta0; p = Y.rows(); dim = p * p; }
-		//! Sample matrix argument
-		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
+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:
+	//! Set internal structures
+	void set_parameters ( const mat &Y0, const double delta0 ) {
+		Y = chmat ( Y0 );
+		delta = delta0;
+		p = Y.rows();
+		dim = p * p;
+	}
+	//! Set internal structures
+	void set_parameters ( const chmat &Y0, const double delta0 ) {
+		Y = Y0;
+		delta = delta0;
+		p = Y.rows();
+		dim = p * p;
+	}
+	//! Sample matrix argument
+	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());
+			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();
 			}
-			//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;}
+		}
+		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;
+	}
 };
 
@@ -1216,91 +1331,95 @@
 /*! Being computed by conversion from `standard' Wishart
 */
-class eiWishartCh: public epdf
-{
-	protected:
-		//! Internal instance of Wishart density
-		eWishartCh W;
-		//! size of Ch
-		int p;
-		//! parameter delta
-		double delta;
-	public:
-		//! constructor function
-		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);}
-		//! access function
-		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
+class eiWishartCh: public epdf {
+protected:
+	//! Internal instance of Wishart density
+	eWishartCh W;
+	//! size of Ch
+	int p;
+	//! parameter delta
+	double delta;
+public:
+	//! constructor function
+	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 );
+	}
+	//! access function
+	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;
-		};
+		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;
+	};
 
 };
 
 //! Random Walk on inverse Wishart
-class rwiWishartCh : public pdf_internal<eiWishartCh>
-{
-	protected:
-		//!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;
-		//! power of the reference
-		double l;
-		//! dimension
-		int p;
-
-	public:
-		rwiWishartCh() : sqd (0), l (0), p (0) {}
-		//! constructor function
-		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);
-
-			iepdf.set_parameters (eye (p), delta);
-			dimc = iepdf.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++;
-			}
-
-			iepdf._setY (sqd*z);
-		}
+class rwiWishartCh : public pdf_internal<eiWishartCh> {
+protected:
+	//!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;
+	//! power of the reference
+	double l;
+	//! dimension
+	int p;
+
+public:
+	rwiWishartCh() : sqd ( 0 ), l ( 0 ), p ( 0 ) {}
+	//! constructor function
+	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 );
+
+		iepdf.set_parameters ( eye ( p ), delta );
+		dimc = iepdf.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++;
+		}
+
+		iepdf._setY ( sqd*z );
+	}
 };
 
@@ -1312,91 +1431,117 @@
 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);};
-		//! Set samples
-		void set_parameters (const Array<vec> &Av) {
-			bdm_assert(Av.size()>0,"Empty samples"); 
-			n = Av.size(); 
-			epdf::set_parameters(Av(0).length());
-			w=1/n*ones(n);
-			samples=Av;
-		};
-		//! 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 vec& _sample(int i) const {return samples(i);};
-		//! 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.
-		//! The vector with indeces of new samples is returned in variable \c index. 
-		void resample ( ivec &index, RESAMPLING_METHOD method = SYSTEMATIC);
-
-		//! Resampling without returning index of new particles. 
-		void resample (RESAMPLING_METHOD method = SYSTEMATIC){ivec ind; resample(ind,method);};
-		
-		//! inherited operation : NOT implemented
-		vec sample() const {
-			bdm_error ("Not implemented");
-			return vec();
-		}
-
-		//! inherited operation : NOT implemented
-		double evallog (const vec &val) const {
-			bdm_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);}
+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 );
+	};
+	//! Set samples
+	void set_parameters ( const Array<vec> &Av ) {
+		bdm_assert ( Av.size() > 0, "Empty samples" );
+		n = Av.size();
+		epdf::set_parameters ( Av ( 0 ).length() );
+		w = 1 / n * ones ( n );
+		samples = Av;
+	};
+	//! 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 vec& _sample ( int i ) const {
+		return samples ( i );
+	};
+	//! 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.
+	//! The vector with indeces of new samples is returned in variable \c index.
+	void resample ( ivec &index, RESAMPLING_METHOD method = SYSTEMATIC );
+
+	//! Resampling without returning index of new particles.
+	void resample ( RESAMPLING_METHOD method = SYSTEMATIC ) {
+		ivec ind;
+		resample ( ind, method );
+	};
+
+	//! inherited operation : NOT implemented
+	vec sample() const {
+		bdm_error ( "Not implemented" );
+		return vec();
+	}
+
+	//! inherited operation : NOT implemented
+	double evallog ( const vec &val ) const {
+		bdm_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 );
 				}
 			}
 		}
+	}
 };
 
@@ -1405,6 +1550,5 @@
 
 template<class sq_T>
-void enorm<sq_T>::set_parameters (const vec &mu0, const sq_T &R0)
-{
+void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) {
 //Fixme test dimensions of mu0 and R0;
 	mu = mu0;
@@ -1414,11 +1558,10 @@
 
 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);
+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);
+	UI::get ( Rtmp, set, "R", UI::compulsory );
 	R = Rtmp; // conversion
 	validate();
@@ -1426,6 +1569,5 @@
 
 template<class sq_T>
-void enorm<sq_T>::dupdate (mat &v, double nu)
-{
+void enorm<sq_T>::dupdate ( mat &v, double nu ) {
 	//
 };
@@ -1437,10 +1579,9 @@
 
 template<class sq_T>
-vec enorm<sq_T>::sample() const
-{
-	vec x (dim);
+vec enorm<sq_T>::sample() const {
+	vec x ( dim );
 #pragma omp critical
-	NorRNG.sample_vector (dim, x);
-	vec smp = R.sqrt_mult (x);
+	NorRNG.sample_vector ( dim, x );
+	vec smp = R.sqrt_mult ( x );
 
 	smp += mu;
@@ -1457,16 +1598,14 @@
 
 template<class sq_T>
-double enorm<sq_T>::evallog_nn (const vec &val) const
-{
+double enorm<sq_T>::evallog_nn ( const vec &val ) const {
 	// 1.83787706640935 = log(2pi)
-	double tmp = -0.5 * (R.invqform (mu - val));// - lognc();
+	double tmp = -0.5 * ( R.invqform ( mu - val ) );// - lognc();
 	return  tmp;
 };
 
 template<class sq_T>
-inline double enorm<sq_T>::lognc () const
-{
+inline double enorm<sq_T>::lognc () const {
 	// 1.83787706640935 = log(2pi)
-	double tmp = 0.5 * (R.cols() * 1.83787706640935 + R.logdet());
+	double tmp = 0.5 * ( R.cols() * 1.83787706640935 + R.logdet() );
 	return tmp;
 };
@@ -1500,8 +1639,7 @@
 
 template<class sq_T>
-shared_ptr<epdf> enorm<sq_T>::marginal ( const RV &rvn ) const
-{
+shared_ptr<epdf> enorm<sq_T>::marginal ( const RV &rvn ) const {
 	enorm<sq_T> *tmp = new enorm<sq_T> ();
-	shared_ptr<epdf> narrow(tmp);
+	shared_ptr<epdf> narrow ( tmp );
 	marginal ( rvn, *tmp );
 	return narrow;
@@ -1509,20 +1647,18 @@
 
 template<class sq_T>
-void enorm<sq_T>::marginal ( const RV &rvn, enorm<sq_T> &target ) const
-{
-	bdm_assert (isnamed(), "rv description is not assigned");
-	ivec irvn = rvn.dataind (rv);
-
-	sq_T Rn (R, irvn);  // select rows and columns of R
+void enorm<sq_T>::marginal ( const RV &rvn, enorm<sq_T> &target ) const {
+	bdm_assert ( isnamed(), "rv description is not assigned" );
+	ivec irvn = rvn.dataind ( rv );
+
+	sq_T Rn ( R, irvn );  // select rows and columns of R
 
 	target.set_rv ( rvn );
-	target.set_parameters (mu (irvn), Rn);
+	target.set_parameters ( mu ( irvn ), Rn );
 }
 
 template<class sq_T>
-shared_ptr<pdf> enorm<sq_T>::condition ( const RV &rvn ) const
-{
+shared_ptr<pdf> enorm<sq_T>::condition ( const RV &rvn ) const {
 	mlnorm<sq_T> *tmp = new mlnorm<sq_T> ();
-	shared_ptr<pdf> narrow(tmp);
+	shared_ptr<pdf> narrow ( tmp );
 	condition ( rvn, *tmp );
 	return narrow;
@@ -1530,18 +1666,17 @@
 
 template<class sq_T>
-void enorm<sq_T>::condition ( const RV &rvn, pdf &target ) const
-{
+void enorm<sq_T>::condition ( const RV &rvn, pdf &target ) const {
 	typedef mlnorm<sq_T> TMlnorm;
 
-	bdm_assert (isnamed(), "rvs are not assigned");
-	TMlnorm &uptarget = dynamic_cast<TMlnorm &>(target);
-
-	RV rvc = rv.subt (rvn);
-	bdm_assert ( (rvc._dsize() + rvn._dsize() == rv._dsize()), "wrong rvn");
+	bdm_assert ( isnamed(), "rvs are not assigned" );
+	TMlnorm &uptarget = dynamic_cast<TMlnorm &> ( target );
+
+	RV rvc = rv.subt ( rvn );
+	bdm_assert ( ( 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);
+	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?
@@ -1550,16 +1685,16 @@
 	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());
-
-	uptarget.set_rv (rvn);
-	uptarget.set_rvc (rvc);
-	uptarget.set_parameters (A, mu1 - A*mu2, R_n);
+	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() );
+
+	uptarget.set_rv ( rvn );
+	uptarget.set_rvc ( rvc );
+	uptarget.set_parameters ( A, mu1 - A*mu2, R_n );
 }
 
@@ -1567,16 +1702,17 @@
 ///////
 template<class sq_T>
-void mgnorm<sq_T >::set_parameters (const shared_ptr<fnc> &g0, const sq_T &R0) {
+void mgnorm<sq_T >::set_parameters ( const shared_ptr<fnc> &g0, const sq_T &R0 ) {
 	g = g0;
-	this->iepdf.set_parameters (zeros (g->dimension()), R0);
+	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);};
+void mgnorm<sq_T >::condition ( const vec &cond ) {
+	this->iepdf._mu() = g->eval ( cond );
+};
 
 //! \todo unify this stuff with to_string()
 template<class sq_T>
-std::ostream &operator<< (std::ostream &os,  mlnorm<sq_T> &ml)
-{
+std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_T> &ml ) {
 	os << "A:" << ml.A << endl;
 	os << "mu:" << ml.mu_const << endl;
Index: library/bdm/stat/merger.h
===================================================================
--- library/bdm/stat/merger.h (revision 693)
+++ library/bdm/stat/merger.h (revision 737)
@@ -86,5 +86,5 @@
 
 	//! Default constructor
-	merger_base () : Npoints(0), Nsources(0), DBG(false), dbg_file(0) {
+	merger_base () : Npoints ( 0 ), Nsources ( 0 ), DBG ( false ), dbg_file ( 0 ) {
 	}
 
@@ -133,5 +133,5 @@
 	}
 	//! Set support points from rectangular grid
-	void set_support ( rectangular_support &Sup) {
+	void set_support ( rectangular_support &Sup ) {
 		Npoints = Sup.points();
 		eSmp.set_parameters ( Npoints, false );
@@ -139,13 +139,13 @@
 		eSmp._w() = ones ( Npoints ) / Npoints; //unifrom size of bins
 		//set samples
-		samples(0)=Sup.first_vec();
-		for (int j=1; j < Npoints; j++ ) {
+		samples ( 0 ) = Sup.first_vec();
+		for ( int j = 1; j < Npoints; j++ ) {
 			samples ( j ) = Sup.next_vec();
 		}
 	}
 	//! Set support points from dicrete grid
-	void set_support ( discrete_support &Sup) {
+	void set_support ( discrete_support &Sup ) {
 		Npoints = Sup.points();
-		eSmp.set_parameters(Sup._Spoints());
+		eSmp.set_parameters ( Sup._Spoints() );
 	}
 	//! set debug file
@@ -311,8 +311,8 @@
 	//!\name Constructors
 	//!@{
-	merger_mix ():Ncoms(0), effss_coef(0), stop_niter(0) { }
-
-	merger_mix ( const Array<shared_ptr<pdf> > &S ):
-		Ncoms(0), effss_coef(0), stop_niter(0) {
+	merger_mix () : Ncoms ( 0 ), effss_coef ( 0 ), stop_niter ( 0 ) { }
+
+	merger_mix ( const Array<shared_ptr<pdf> > &S ) :
+			Ncoms ( 0 ), effss_coef ( 0 ), stop_niter ( 0 ) {
 		set_sources ( S );
 	}
Index: library/bdm/stat/discrete.h
===================================================================
--- library/bdm/stat/discrete.h (revision 730)
+++ library/bdm/stat/discrete.h (revision 737)
@@ -62,8 +62,10 @@
 	//! Get next active vector, call ONLY after first_vector()!
 	const vec& next_vec();
-	
+
 	//! return active vector, call ONLY after first_vector()!
-	const vec& act_vec(){return actvec;};
-	
+	const vec& act_vec() {
+		return actvec;
+	};
+
 	//! \todo to je asi navic .. v predkovi!
 	ivec nearest_point ( const vec &val );
@@ -143,8 +145,8 @@
 	//! constructor function fills values by calling function \c f , double f(vec&), given by a pointer
 	void set_values ( double ( *evalptr ) ( const vec& ) );
-	
+
 	//! constructor function fills values by calling function \c f , double f(vec&), given by a pointer
-	void set_values ( const epdf &ep);
-	
+	void set_values ( const epdf &ep );
+
 	//! get value nearest to the given point
 	double nearest_val ( const vec &val ) {
@@ -155,5 +157,7 @@
 		return vec_1 ( nearest_val ( val ) );
 	}
-	const vec & _values() const { return values; }
+	const vec & _values() const {
+		return values;
+	}
 };
 UIREGISTER ( grid_fnc );
Index: library/bdm/stat/emix.h
===================================================================
--- library/bdm/stat/emix.h (revision 716)
+++ library/bdm/stat/emix.h (revision 737)
@@ -53,10 +53,10 @@
 		// adjust rv and rvc
 
-		set_rv( rv ); // TODO co kdyby tohle samo uz nastavovalo dimension?!?!
+		set_rv ( rv ); // TODO co kdyby tohle samo uz nastavovalo dimension?!?!
 		dim = rv._dsize();
 
 		rvc = nom0->_rv().subt ( rv );
 		dimc = rvc._dsize();
-	
+
 		//prepare data structures
 		if ( copy ) {
@@ -95,5 +95,5 @@
 	// not implemented
 	mratio ( const mratio & );
-	mratio &operator=( const mratio & );
+	mratio &operator= ( const mratio & );
 };
 
@@ -120,10 +120,10 @@
 	emix ( ) : epdf ( ) { }
 
-	  /*!
-	    \brief Set weights \c w and components \c Coms
-
-	    Shared pointers in Coms are kept inside this instance and
-	    shouldn't be modified after being passed to this method.
-	  */
+	/*!
+	  \brief Set weights \c w and components \c Coms
+
+	  Shared pointers in Coms are kept inside this instance and
+	  shouldn't be modified after being passed to this method.
+	*/
 	void set_parameters ( const vec &w, const Array<shared_ptr<epdf> > &Coms );
 
@@ -209,22 +209,21 @@
 	//!@}
 	void from_setting ( const Setting &set ) {
-	
-		vec w0;			
+
+		vec w0;
 		Array<shared_ptr<epdf> > Coms0;
-	
+
 		UI::get ( Coms0, set, "pdfs", UI::compulsory );
 
-		if( !UI::get( w0, set, "weights", UI::optional ) )
-		{
+		if ( !UI::get ( w0, set, "weights", UI::optional ) ) {
 			int len = Coms.length();
-			w0.set_length( len );
+			w0.set_length ( len );
 			w0 = 1.0 / len;
 		}
 
 		// TODO asi lze nacitat primocare do w a coms, jen co bude hotovy validate()
-		set_parameters( w0, Coms0 );
+		set_parameters ( w0, Coms0 );
 	}
 };
-SHAREDPTR( emix );
+SHAREDPTR ( emix );
 UIREGISTER ( emix );
 
@@ -333,5 +332,5 @@
 	}
 	//! Set internal \c pdfs from given values
-	void set_elements (const Array<shared_ptr<pdf> > &mFacs );
+	void set_elements ( const Array<shared_ptr<pdf> > &mFacs );
 
 	double evallogcond ( const vec &val, const vec &cond ) {
@@ -375,5 +374,5 @@
 		for ( int i = ( pdfs.length() - 1 ); i >= 0; i-- ) {
 			// generate contribution of this pdf
-			smpi = pdfs(i)->samplecond(dls ( i )->get_cond ( smp , cond ));			
+			smpi = pdfs ( i )->samplecond ( dls ( i )->get_cond ( smp , cond ) );
 			// copy contribution of this pdf into smp
 			dls ( i )->pushup ( smp, smpi );
@@ -408,5 +407,5 @@
 	//! Default constructor
 	eprod () : epdfs ( 0 ), dls ( 0 ) {};
-	//! Set internal 
+	//! Set internal
 	void set_parameters ( const Array<const epdf*> &epdfs0, bool named = true ) {
 		epdfs = epdfs0;//.set_length ( epdfs0.length() );
@@ -500,5 +499,5 @@
 public:
 	//!Default constructor
-	mmix() : Coms(0) { }
+	mmix() : Coms ( 0 ) { }
 
 	//! Set weights \c w and components \c R
@@ -506,22 +505,22 @@
 		//!\todo check if all components are OK
 		Coms = Coms0;
-		w=w0;	
-
-		if (Coms0.length()>0){
-			set_rv(Coms(0)->_rv());
+		w = w0;
+
+		if ( Coms0.length() > 0 ) {
+			set_rv ( Coms ( 0 )->_rv() );
 			dim = rv._dsize();
-			set_rvc(Coms(0)->_rvc());
+			set_rvc ( Coms ( 0 )->_rvc() );
 			dimc = rvc._dsize();
 		}
 	}
-	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);
+	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);
+	vec samplecond ( const vec &cond );
 
 	//! Load from structure with elements:
@@ -534,23 +533,22 @@
 	//!@}
 	void from_setting ( const Setting &set ) {
-		UI::get ( Coms, set, "pdfs", UI::compulsory );	
+		UI::get ( Coms, set, "pdfs", UI::compulsory );
 
 		// TODO ma byt zde, ci ve validate()?
-		if (Coms.length()>0){
-			set_rv(Coms(0)->_rv());
+		if ( Coms.length() > 0 ) {
+			set_rv ( Coms ( 0 )->_rv() );
 			dim = rv._dsize();
-			set_rvc(Coms(0)->_rvc());
+			set_rvc ( Coms ( 0 )->_rvc() );
 			dimc = rvc._dsize();
 		}
 
-		if( !UI::get( w, set, "weights", UI::optional ) )
-		{
+		if ( !UI::get ( w, set, "weights", UI::optional ) ) {
 			int len = Coms.length();
-			w.set_length( len );
+			w.set_length ( len );
 			w = 1.0 / len;
 		}
 	}
 };
-SHAREDPTR( mmix );
+SHAREDPTR ( mmix );
 UIREGISTER ( mmix );
 
Index: library/bdm/stat/exp_family.cpp
===================================================================
--- library/bdm/stat/exp_family.cpp (revision 730)
+++ library/bdm/stat/exp_family.cpp (revision 737)
@@ -12,17 +12,17 @@
 using std::cout;
 
-	///////////
-
-void BMEF::bayes( const vec &yt, const vec &cond ) {
+///////////
+
+void BMEF::bayes ( const vec &yt, const vec &cond ) {
 	this->bayes_weighted ( yt, cond, 1.0 );
 };
 
-void egiw::set_parameters (int dimx0, ldmat V0, double nu0) {
+void egiw::set_parameters ( int dimx0, ldmat V0, double nu0 ) {
 	dimx = dimx0;
 	nPsi = V0.rows() - dimx;
-	dim = dimx * (dimx + nPsi); // size(R) + size(Theta)
-	
+	dim = dimx * ( dimx + nPsi ); // size(R) + size(Theta)
+
 	V = V0;
-	if (nu0 < 0) {
+	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
@@ -35,62 +35,62 @@
 	mat M;
 	chmat R;
-	sample_mat(M,R);
-	
-	return concat (cvectorize(M),cvectorize(R.to_mat()));
-}
-
-mat egiw::sample_mat(int n) const {
+	sample_mat ( M, R );
+
+	return concat ( cvectorize ( M ), cvectorize ( R.to_mat() ) );
+}
+
+mat egiw::sample_mat ( int n ) const {
 	// TODO - correct approach - convert to product of norm * Wishart
 	mat M;
 	ldmat Vz;
 	ldmat Lam;
-	factorize(M,Vz,Lam);
-	
-	chmat ChLam(Lam.to_mat());
+	factorize ( M, Vz, Lam );
+
+	chmat ChLam ( Lam.to_mat() );
 	chmat iChLam;
-	ChLam.inv(iChLam);
-	
-	eWishartCh Omega; //inverse Wishart, result is R, 
-	Omega.set_parameters(iChLam,nu-2*nPsi-dimx); // 2*nPsi is there to match numercial simulations - check if analytically correct
-	
+	ChLam.inv ( iChLam );
+
+	eWishartCh Omega; //inverse Wishart, result is R,
+	Omega.set_parameters ( iChLam, nu - 2*nPsi - dimx ); // 2*nPsi is there to match numercial simulations - check if analytically correct
+
 	mat OmChi;
-	mat Z(M.rows(),M.cols());
+	mat Z ( M.rows(), M.cols() );
 
 	mat Mi;
 	mat RChiT;
-	mat tmp(dimension(), n);
-	for (int i=0; i<n;i++){
-		OmChi=Omega.sample_mat();
-		RChiT=inv(OmChi);
-		Z=randn(M.rows(), M.cols());
-		Mi = M + RChiT * Z * inv(Vz._L().T() *diag(sqrt(Vz._D())));
-	
-		tmp.set_col(i,concat (cvectorize(Mi),cvectorize(RChiT*RChiT.T())));
+	mat tmp ( dimension(), n );
+	for ( int i = 0; i < n; i++ ) {
+		OmChi = Omega.sample_mat();
+		RChiT = inv ( OmChi );
+		Z = randn ( M.rows(), M.cols() );
+		Mi = M + RChiT * Z * inv ( Vz._L().T() * diag ( sqrt ( Vz._D() ) ) );
+
+		tmp.set_col ( i, concat ( cvectorize ( Mi ), cvectorize ( RChiT*RChiT.T() ) ) );
 	}
 	return tmp;
 }
 
-void egiw::sample_mat(mat &Mi, chmat &Ri)const{
-	
+void egiw::sample_mat ( mat &Mi, chmat &Ri ) const {
+
 	// TODO - correct approach - convert to product of norm * Wishart
 	mat M;
 	ldmat Vz;
 	ldmat Lam;
-	factorize(M,Vz,Lam);
-	
+	factorize ( M, Vz, Lam );
+
 	chmat Ch;
-	Ch.setCh(Lam._L()*diag(sqrt(Lam._D())));
+	Ch.setCh ( Lam._L() *diag ( sqrt ( Lam._D() ) ) );
 	chmat iCh;
-	Ch.inv(iCh);
-	
-	eWishartCh Omega; //inverse Wishart, result is R, 
-	Omega.set_parameters(iCh,nu-2*nPsi-dimx); // 2*nPsi is there to match numercial simulations - check if analytically correct
-	
+	Ch.inv ( iCh );
+
+	eWishartCh Omega; //inverse Wishart, result is R,
+	Omega.set_parameters ( iCh, nu - 2*nPsi - dimx ); // 2*nPsi is there to match numercial simulations - check if analytically correct
+
 	chmat Omi;
-	Omi.setCh(Omega.sample_mat());
-	
-	mat Z=randn(M.rows(), M.cols());
-	Mi = M+ Omi._Ch() * Z * inv(Vz._L()*diag(sqrt(Vz._D())));
-	Omi.inv(Ri);
+	Omi.setCh ( Omega.sample_mat() );
+
+	mat Z = randn ( M.rows(), M.cols() );
+	Mi = M + Omi._Ch() * Z * inv ( Vz._L() * diag ( sqrt ( Vz._D() ) ) );
+	Omi.inv ( Ri );
 }
 
@@ -140,6 +140,6 @@
 
 //	bdm_assert_debug ( ( ( -nkG - nkW ) > -Inf ) && ( ( -nkG - nkW ) < Inf ), "ARX improper" );
-	if (-nkG - nkW==Inf){
-		cout << "??" <<endl;
+	if ( -nkG - nkW == Inf ) {
+		cout << "??" << endl;
 	}
 	return -nkG - nkW;
@@ -162,24 +162,24 @@
 }
 
-void egiw::factorize(mat &M, ldmat &Vz, ldmat &Lam) const{	
+void egiw::factorize ( mat &M, ldmat &Vz, ldmat &Lam ) const {
 	const mat &L = V._L();
 	const vec &D = V._D();
 	int end = L.rows() - 1;
-	
-	Vz=ldmat(L ( dimx, end, dimx, end ), D(dimx,end));
-	mat iLsub = ltuinv ( Vz._L());
+
+	Vz = ldmat ( L ( dimx, end, dimx, end ), D ( dimx, end ) );
+	mat iLsub = ltuinv ( Vz._L() );
 	// set mean value
 	mat Lpsi = L ( dimx, end, 0, dimx - 1 );
 	M = iLsub * Lpsi;
-	
-	Lam =ldmat ( L (0, dimx-1, 0, dimx-1 ), D (0, dimx-1 ) );  //exp val of R
-	 if (1){ // test with Peterka
-		 mat VF=V.to_mat();
-		 mat Vf=VF(0,dimx-1,0, dimx-1);
-		 mat Vzf = VF(dimx,end,0,dimx-1);
-		 mat VZ = VF(dimx,end,dimx,end);
-		 
-		 mat Lam2 = Vf-Vzf.T()*inv(VZ)*Vzf;
-	 }
+
+	Lam = ldmat ( L ( 0, dimx - 1, 0, dimx - 1 ), D ( 0, dimx - 1 ) );  //exp val of R
+	if ( 1 ) { // test with Peterka
+		mat VF = V.to_mat();
+		mat Vf = VF ( 0, dimx - 1, 0, dimx - 1 );
+		mat Vzf = VF ( dimx, end, 0, dimx - 1 );
+		mat VZ = VF ( dimx, end, dimx, end );
+
+		mat Lam2 = Vf - Vzf.T() * inv ( VZ ) * Vzf;
+	}
 }
 
@@ -193,5 +193,5 @@
 //		mat Dsub = diag ( D ( 1, end ) );
 
-		ldmat LD(inv(Lsub).T(), 1.0/D(1,end));
+		ldmat LD ( inv ( Lsub ).T(), 1.0 / D ( 1, end ) );
 		return LD;
 
@@ -217,5 +217,5 @@
 		mat R;
 		mean_mat ( M, R );
-		return concat( cvectorize (  M),cvectorize( R ) );
+		return concat ( cvectorize ( M ), cvectorize ( R ) );
 	}
 
@@ -229,10 +229,10 @@
 	ldmat itmp ( l );
 	tmp.inv ( itmp );
-	
+
 	// following Wikipedia notation
-	// m=nu-nPsi-dimx-1, p=dimx 
-	double mp1p=nu-nPsi-2*dimx; // m-p+1
-	double mp1m=mp1p-2;         // m-p-1
-	
+	// m=nu-nPsi-dimx-1, p=dimx
+	double mp1p = nu - nPsi - 2 * dimx; // m-p+1
+	double mp1m = mp1p - 2;     // m-p-1
+
 	if ( dimx == 1 ) {
 		double cove = V._D() ( 0 ) / mp1m ;
@@ -240,28 +240,28 @@
 		vec var ( l );
 		var.set_subvector ( 0, diag ( itmp.to_mat() ) *cove );
-		var ( l - 1 ) = cove * cove / (mp1m-2);
+		var ( l - 1 ) = cove * cove / ( mp1m - 2 );
 		return var;
 	} else {
-		ldmat Vll( V, linspace(0,dimx-1)); // top-left part of V
-		mat Y=Vll.to_mat();
-		mat varY(Y.rows(), Y.cols()); 
-		
-		double denom = (mp1p-1)*mp1m*mp1m*(mp1m-2);         // (m-p)(m-p-1)^2(m-p-3)
-		
-		int i,j;
-		for ( i=0; i<Y.rows(); i++){
-			for ( j=0; j<Y.cols(); j++){
-				varY(i,j) = (mp1p*Y(i,j)*Y(i,j) + mp1m * Y(i,i)* Y(j,j)) /denom;
+		ldmat Vll ( V, linspace ( 0, dimx - 1 ) ); // top-left part of V
+		mat Y = Vll.to_mat();
+		mat varY ( Y.rows(), Y.cols() );
+
+		double denom = ( mp1p - 1 ) * mp1m * mp1m * ( mp1m - 2 );         // (m-p)(m-p-1)^2(m-p-3)
+
+		int i, j;
+		for ( i = 0; i < Y.rows(); i++ ) {
+			for ( j = 0; j < Y.cols(); j++ ) {
+				varY ( i, j ) = ( mp1p * Y ( i, j ) * Y ( i, j ) + mp1m * Y ( i, i ) * Y ( j, j ) ) / denom;
 			}
 		}
-		vec mean_dR = diag(Y)/mp1m; // corresponds to cove
-		vec var_th=diag ( itmp.to_mat() );
-		vec var_Th ( mean_dR.length()*var_th.length() );
+		vec mean_dR = diag ( Y ) / mp1m; // corresponds to cove
+		vec var_th = diag ( itmp.to_mat() );
+		vec var_Th ( mean_dR.length() *var_th.length() );
 		// diagonal of diag(mean_dR) \kron diag(var_th)
-		for (int i=0; i<mean_dR.length(); i++){
-			var_Th.set_subvector(i*var_th.length(), var_th*mean_dR(i));  
-		}
-		
-		return concat(var_Th, cvectorize(varY));
+		for ( int i = 0; i < mean_dR.length(); i++ ) {
+			var_Th.set_subvector ( i*var_th.length(), var_th*mean_dR ( i ) );
+		}
+
+		return concat ( var_Th, cvectorize ( varY ) );
 	}
 }
