Index: library/bdm/stat/exp_family.h
===================================================================
--- library/bdm/stat/exp_family.h (revision 1063)
+++ library/bdm/stat/exp_family.h (revision 1064)
@@ -37,37 +37,37 @@
 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 NOT_IMPLEMENTED(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 ) NOT_IMPLEMENTED_VOID;
+    //! 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 NOT_IMPLEMENTED(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 ) NOT_IMPLEMENTED_VOID;
 };
 
@@ -75,48 +75,50 @@
 //! Estimator for Exponential family
 class BMEF : public BM {
-	public:
-	//! forgetting factor
-	double frg;
-	protected:
-	//! cached value of lognc() in the previous step (used in evaluation of \c ll )
-	double last_lognc;
-	//! factor k = [0..1] for scheduling of forgetting factor: \f$ frg_t = (1-k) * frg_{t-1} + k \f$, default 0
-	double frg_sched_factor;
-public:
-	//! Default constructor (=empty constructor)
-	BMEF ( double frg0 = 1.0 ) : BM (), frg ( frg0 ), last_lognc(0.0),frg_sched_factor(0.0) {}
-	//! Copy constructor
-	BMEF ( const BMEF &B ) : BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ),frg_sched_factor(B.frg_sched_factor) {}
-	//!get statistics from another model
-	virtual void set_statistics ( const BMEF* BM0 ) NOT_IMPLEMENTED_VOID;
-
-	//! Weighted update of sufficient statistics (Bayes rule)
-	virtual void bayes_weighted ( const vec &data, const vec &cond = empty_vec, const double w = 1.0 ) {
-		if (frg_sched_factor>0) {frg = frg*(1-frg_sched_factor)+frg_sched_factor;}
-	};
-	//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, double weight=1.0 ) NOT_IMPLEMENTED_VOID;;
-
-
-	void to_setting ( Setting &set ) const
-	{			
-		BM::to_setting( set );
-		UI::save(frg, set, "frg");
-		UI::save( frg_sched_factor, set, "frg_sched_factor" );
-	} 
-
-	void from_setting( const Setting &set) {
-		BM::from_setting(set);
-		if ( !UI::get ( frg, set, "frg" ) )
-			frg = 1.0;
-		UI::get ( frg_sched_factor, set, "frg_sched_factor",UI::optional );
-	}
-
-	void validate() {
-		BM::validate();	
-	}
+public:
+    //! forgetting factor
+    double frg;
+protected:
+    //! cached value of lognc() in the previous step (used in evaluation of \c ll )
+    double last_lognc;
+    //! factor k = [0..1] for scheduling of forgetting factor: \f$ frg_t = (1-k) * frg_{t-1} + k \f$, default 0
+    double frg_sched_factor;
+public:
+    //! Default constructor (=empty constructor)
+    BMEF ( double frg0 = 1.0 ) : BM (), frg ( frg0 ), last_lognc(0.0),frg_sched_factor(0.0) {}
+    //! Copy constructor
+    BMEF ( const BMEF &B ) : BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ),frg_sched_factor(B.frg_sched_factor) {}
+    //!get statistics from another model
+    virtual void set_statistics ( const BMEF* BM0 ) NOT_IMPLEMENTED_VOID;
+
+    //! Weighted update of sufficient statistics (Bayes rule)
+    virtual void bayes_weighted ( const vec &data, const vec &cond = empty_vec, const double w = 1.0 ) {
+        if (frg_sched_factor>0) {
+            frg = frg*(1-frg_sched_factor)+frg_sched_factor;
+        }
+    };
+    //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, double weight=1.0 ) NOT_IMPLEMENTED_VOID;;
+
+
+    void to_setting ( Setting &set ) const
+    {
+        BM::to_setting( set );
+        UI::save(frg, set, "frg");
+        UI::save( frg_sched_factor, set, "frg_sched_factor" );
+    }
+
+    void from_setting( const Setting &set) {
+        BM::from_setting(set);
+        if ( !UI::get ( frg, set, "frg" ) )
+            frg = 1.0;
+        UI::get ( frg_sched_factor, set, "frg_sched_factor",UI::optional );
+    }
+
+    void validate() {
+        BM::validate();
+    }
 
 };
@@ -127,19 +129,19 @@
 where \f$ x_t \f$ is the \c rv, \f$ y_t \f$ is the \c rvc and g is a deterministic transformation of class fn.
 */
-class mgdirac: public pdf{
-	protected:
-	shared_ptr<fnc> g;
-	public:
-		vec samplecond(const vec &cond) {
-			bdm_assert_debug(cond.length()==g->dimensionc(),"given cond in not compatible with g");
-			vec tmp = g->eval(cond);
-			return tmp;
-		} 
-		double evallogcond ( const vec &yt, const vec &cond ){
-			return std::numeric_limits< double >::max();
-		}
-		void from_setting(const Setting& set);
-		void to_setting(Setting &set) const;
-		void validate();
+class mgdirac: public pdf {
+protected:
+    shared_ptr<fnc> g;
+public:
+    vec samplecond(const vec &cond) {
+        bdm_assert_debug(cond.length()==g->dimensionc(),"given cond in not compatible with g");
+        vec tmp = g->eval(cond);
+        return tmp;
+    }
+    double evallogcond ( const vec &yt, const vec &cond ) {
+        return std::numeric_limits< double >::max();
+    }
+    void from_setting(const Setting& set);
+    void to_setting(Setting &set) const;
+    void validate();
 };
 UIREGISTER(mgdirac);
@@ -157,93 +159,93 @@
 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 to_setting ( Setting &root ) const ;
-	
-	void validate();
-	//!@}
-
-	//! \name Mathematical operations
-	//!@{
-
-	//! dupdate in exponential form (not really handy)
-	void dupdate ( mat &v, double nu = 1.0 );
-
-	//! evaluate bhattacharya distance
-	double bhattacharyya(const enorm<sq_T> &e2){
-		bdm_assert(dim == e2.dimension(), "enorms of differnt dimensions");
-		sq_T P=R;
-		P.add(e2._R());
-		
-		double tmp = 0.125*P.invqform(mu - e2._mu()) + 0.5*(P.logdet() - 0.5*(R.logdet() + e2._R().logdet()));
-		return tmp;
-	}
-	
-	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() );
-	}
-	mat covariance() const {
-		return 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;
-	}
-	//!@}
+    //! 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 to_setting ( Setting &root ) const ;
+
+    void validate();
+    //!@}
+
+    //! \name Mathematical operations
+    //!@{
+
+    //! dupdate in exponential form (not really handy)
+    void dupdate ( mat &v, double nu = 1.0 );
+
+    //! evaluate bhattacharya distance
+    double bhattacharyya(const enorm<sq_T> &e2) {
+        bdm_assert(dim == e2.dimension(), "enorms of differnt dimensions");
+        sq_T P=R;
+        P.add(e2._R());
+
+        double tmp = 0.125*P.invqform(mu - e2._mu()) + 0.5*(P.logdet() - 0.5*(R.logdet() + e2._R().logdet()));
+        return tmp;
+    }
+
+    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() );
+    }
+    mat covariance() const {
+        return 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;
+    }
+    //!@}
 
 };
@@ -255,5 +257,5 @@
 SHAREDPTR2 ( enorm, fsqmat );
 
-//! \class bdm::egauss 
+//! \class bdm::egauss
 //!\brief Gaussian (Normal) distribution. Same as enorm<fsqmat>.
 typedef enorm<ldmat> egauss;
@@ -266,72 +268,82 @@
 /*! distribution of multivariate Student t density
 
-Based on article by Genest and Zidek, 
+Based on article by Genest and Zidek,
 */
 template<class sq_T>
-class estudent : public eEF{
-	protected:
-		//! mena value
-		vec mu;
-		//! matrix H
-		sq_T H;
-		//! degrees of freedom
-		double delta;
-	public:
-		double evallog_nn(const vec &val) const{
-			double tmp = -0.5*H.logdet() - 0.5*(delta + dim) * log(1+ H.invqform(val - mu)/delta);
-			return tmp;
-		}
-		double lognc() const {
-			//log(pi) = 1.14472988584940
-			double tmp = -lgamma(0.5*(delta+dim))+lgamma(0.5*delta) + 0.5*dim*(log(delta) + 1.14472988584940);
-			return tmp;
-		}
-		void marginal (const RV &rvm, estudent<sq_T> &marg) const {
-			ivec ind = rvm.findself_ids(rv); // indices of rvm in rv
-			marg._mu() = mu(ind);
-			marg._H() = sq_T(H,ind);
-			marg._delta() = delta;
-			marg.validate();
-		}
-		shared_ptr<epdf> marginal(const RV &rvm) const {
-			shared_ptr<estudent<sq_T> > tmp = new estudent<sq_T>;
-			marginal(rvm, *tmp);
-			return tmp;
-		}
-		vec sample() const NOT_IMPLEMENTED(vec(0))
-		
-		vec mean() const {return mu;}
-		mat covariance() const {
-			return delta/(delta-2)*H.to_mat();
-		}
-		vec variance() const {return diag(covariance());}
-			//! \name access
-		//! @{
-			//! access function
-			vec& _mu() {return mu;}
-			//! access function
-			sq_T& _H() {return H;}
-			//! access function
-			double& _delta() {return delta;}
-			//!@}
-			//! todo
-			void from_setting(const Setting &set){
-				epdf::from_setting(set);
-				mat H0;
-				UI::get(H0,set, "H");
-				H= H0; // conversion!!
-				UI::get(delta,set,"delta");
-				UI::get(mu,set,"mu");
-			}
-			void to_setting(Setting &set) const{
-				epdf::to_setting(set);
-				UI::save(H.to_mat(), set, "H");
-				UI::save(delta, set, "delta");
-				UI::save(mu, set, "mu");
-			}
-			void validate() {
-				eEF::validate();
-				dim = H.rows();
-			}
+class estudent : public eEF {
+protected:
+    //! mena value
+    vec mu;
+    //! matrix H
+    sq_T H;
+    //! degrees of freedom
+    double delta;
+public:
+    double evallog_nn(const vec &val) const {
+        double tmp = -0.5*H.logdet() - 0.5*(delta + dim) * log(1+ H.invqform(val - mu)/delta);
+        return tmp;
+    }
+    double lognc() const {
+        //log(pi) = 1.14472988584940
+        double tmp = -lgamma(0.5*(delta+dim))+lgamma(0.5*delta) + 0.5*dim*(log(delta) + 1.14472988584940);
+        return tmp;
+    }
+    void marginal (const RV &rvm, estudent<sq_T> &marg) const {
+        ivec ind = rvm.findself_ids(rv); // indices of rvm in rv
+        marg._mu() = mu(ind);
+        marg._H() = sq_T(H,ind);
+        marg._delta() = delta;
+        marg.validate();
+    }
+    shared_ptr<epdf> marginal(const RV &rvm) const {
+        shared_ptr<estudent<sq_T> > tmp = new estudent<sq_T>;
+        marginal(rvm, *tmp);
+        return tmp;
+    }
+    vec sample() const NOT_IMPLEMENTED(vec(0))
+
+    vec mean() const {
+        return mu;
+    }
+    mat covariance() const {
+        return delta/(delta-2)*H.to_mat();
+    }
+    vec variance() const {
+        return diag(covariance());
+    }
+    //! \name access
+    //! @{
+    //! access function
+    vec& _mu() {
+        return mu;
+    }
+    //! access function
+    sq_T& _H() {
+        return H;
+    }
+    //! access function
+    double& _delta() {
+        return delta;
+    }
+    //!@}
+    //! todo
+    void from_setting(const Setting &set) {
+        epdf::from_setting(set);
+        mat H0;
+        UI::get(H0,set, "H");
+        H= H0; // conversion!!
+        UI::get(delta,set,"delta");
+        UI::get(mu,set,"mu");
+    }
+    void to_setting(Setting &set) const {
+        epdf::to_setting(set);
+        UI::save(H.to_mat(), set, "H");
+        UI::save(delta, set, "delta");
+        UI::save(mu, set, "mu");
+    }
+    void validate() {
+        eEF::validate();
+        dim = H.rows();
+    }
 };
 UIREGISTER2(estudent,fsqmat);
@@ -346,122 +358,122 @@
 */
 class egiw : public eEF {
-	//! \var log_level_enums logvartheta
-	//! Log variance of the theta part 
-
-	LOG_LEVEL(egiw,logvartheta);
-
-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(),dimx(0) {};
-	egiw ( int dimx0, ldmat V0, double nu0 = -1.0 ) : eEF(),dimx(0) {
-		set_parameters ( dimx0, V0, nu0 );
-		validate();
-	};
-
-	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;
-	//mat covariance() 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;
-	};
-
-	//! marginal density (only student for now)
-	shared_ptr<epdf> marginal(const RV &rvm) const {
-		bdm_assert(dimx==1, "Not supported");
-		//TODO - this is too trivial!!!
-		ivec ind = rvm.findself_ids(rv);
-		if (min(ind)==0) { //assume it si 
-			shared_ptr<estudent<ldmat> > tmp = new estudent<ldmat>;
-			mat M;
-			ldmat Vz;
-			ldmat Lam;
-			factorize(M,Vz,Lam);
-			
-			tmp->_mu() = M.get_col(0);
-			ldmat H;
-			Vz.inv(H);
-			H *=Lam._D()(0)/nu;
-			tmp->_H() = H;
-			tmp->_delta() = nu;
-			tmp->validate();
-			return tmp;
-		}
-		return NULL;
-	}
-	//! \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.L     = [];             // L part of matrix V
-	V.D     = [];             // D part of matrix V
-	-or- V  = []              // full matrix V
-	-or- 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
-
-	\sa log_level_enums
-	*/
-	void from_setting ( const Setting &set );
-	//! see egiw::from_setting
-	void to_setting ( Setting& set ) const;
-	void validate();
-	void log_register ( bdm::logger& L, const string& prefix );
-
-	void log_write() const;
-	//!@}
+    //! \var log_level_enums logvartheta
+    //! Log variance of the theta part
+
+    LOG_LEVEL(egiw,logvartheta);
+
+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(),dimx(0) {};
+    egiw ( int dimx0, ldmat V0, double nu0 = -1.0 ) : eEF(),dimx(0) {
+        set_parameters ( dimx0, V0, nu0 );
+        validate();
+    };
+
+    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;
+    //mat covariance() 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;
+    };
+
+    //! marginal density (only student for now)
+    shared_ptr<epdf> marginal(const RV &rvm) const {
+        bdm_assert(dimx==1, "Not supported");
+        //TODO - this is too trivial!!!
+        ivec ind = rvm.findself_ids(rv);
+        if (min(ind)==0) { //assume it si
+            shared_ptr<estudent<ldmat> > tmp = new estudent<ldmat>;
+            mat M;
+            ldmat Vz;
+            ldmat Lam;
+            factorize(M,Vz,Lam);
+
+            tmp->_mu() = M.get_col(0);
+            ldmat H;
+            Vz.inv(H);
+            H *=Lam._D()(0)/nu;
+            tmp->_H() = H;
+            tmp->_delta() = nu;
+            tmp->validate();
+            return tmp;
+        }
+        return NULL;
+    }
+    //! \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.L     = [];             // L part of matrix V
+    V.D     = [];             // D part of matrix V
+    -or- V  = []              // full matrix V
+    -or- 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
+
+    \sa log_level_enums
+    */
+    void from_setting ( const Setting &set );
+    //! see egiw::from_setting
+    void to_setting ( Setting& set ) const;
+    void validate();
+    void log_register ( bdm::logger& L, const string& prefix );
+
+    void log_write() const;
+    //!@}
 };
 UIREGISTER ( egiw );
@@ -478,79 +490,79 @@
 class eDirich: public eEF {
 protected:
-	//!sufficient statistics
-	vec beta;
-public:
-	//!\name Constructors
-	//!@{
-
-	eDirich () : eEF () {};
-	eDirich ( const eDirich &D0 ) : eEF () {
-		set_parameters ( D0.beta );
-		validate();
-	};
-	eDirich ( const vec &beta0 ) {
-		set_parameters ( beta0 );
-		validate();
-	};
-	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 );
+    //!sufficient statistics
+    vec beta;
+public:
+    //!\name Constructors
+    //!@{
+
+    eDirich () : eEF () {};
+    eDirich ( const eDirich &D0 ) : eEF () {
+        set_parameters ( D0.beta );
+        validate();
+    };
+    eDirich ( const vec &beta0 ) {
+        set_parameters ( beta0 );
+        validate();
+    };
+    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;
-	}
-
-	/*! Create object from the following structure
-	\code
-	class = 'eDirich';
-	beta  = [...];           % vector parameter beta
-	--- inherited fields ---
-	bdm::eEF::from_setting
-	\endcode
-	*/
-	void from_setting ( const Setting &set );
-
-	void validate();
-
-	void to_setting ( Setting &set ) const;
+            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;
+    }
+
+    /*! Create object from the following structure
+    \code
+    class = 'eDirich';
+    beta  = [...];           % vector parameter beta
+    --- inherited fields ---
+    bdm::eEF::from_setting
+    \endcode
+    */
+    void from_setting ( const Setting &set );
+
+    void validate();
+
+    void to_setting ( Setting &set ) const;
 };
 UIREGISTER ( eDirich );
@@ -560,88 +572,88 @@
 Continuous Dirichlet density of \f$n\f$-dimensional variable \f$x\f$
 \f[
-f(x|\alpha,\beta)  = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\, x^{\alpha-1}(1-x)^{\beta-1} 
+f(x|\alpha,\beta)  = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\, x^{\alpha-1}(1-x)^{\beta-1}
 \f]
 is a simplification of Dirichlet to univariate case.
 */
 class eBeta: public eEF {
-	public:
-		//!sufficient statistics
-		vec alpha;
-		//!sufficient statistics
-		vec beta;
-	public:
-	//!\name Constructors
-	//!@{
-			
-	eBeta () : eEF () {};
-	eBeta ( const eBeta &B0 ) : eEF (), alpha(B0.alpha),beta(B0.beta) {
-		validate();
-	};
-	//!@}
-			
-	//! using sampling procedure from wikipedia
-	vec sample() const {
-		vec y ( beta.length() ); // all vectors
-		for ( int i = 0; i < beta.length(); i++ ) {
-			GamRNG.setup ( alpha ( i ), 1 );
-			#pragma omp critical
-			double Ga = GamRNG();
-					
-			GamRNG.setup ( beta ( i ), 1 );
-			#pragma omp critical
-			double Gb = GamRNG();
-					
-			y ( i ) = Ga/(Ga+Gb);
-		}
-		return y;
-	}
-			
-	vec mean() const {
-		return elem_div(alpha, alpha + beta); // dot-division
-	};
-	vec variance() const {
-		vec apb=alpha+beta;
-		return elem_div (elem_mult ( alpha, beta) ,
-							elem_mult ( elem_mult(apb,apb), apb+1 ) );
-	}
-	//! In this instance, val is ...
-	double evallog_nn ( const vec &val ) const {
-		double tmp;
-		tmp = ( alpha - 1 ) * log ( val ) + (beta-1)*log(1-val);
-		return tmp;
-	}
-			
-	double lognc () const {
-		double lgb = 0.0;
-		for ( int i = 0; i < beta.length(); i++ ) {
-			lgb += -lgamma ( alpha(i)+beta(i) ) + lgamma(alpha(i)) + lgamma(beta(i));
-		}
-		return lgb;
-	}
-			
-	/*! Create object from the following structure
-
-	\code
-	class = 'eBeta';
-	alpha = [...];           % vector parameter alpha
-	beta  = [...];           % vector parameter beta of the same length as alpha
-	\endcode
-			
-	Class does not call bdm::eEF::from_setting
-	*/
-	void from_setting ( const Setting &set ){
-		UI::get(alpha, set, "alpha", UI::compulsory);
-		UI::get(beta, set, "beta", UI::compulsory);
-	}
-
-	void validate(){
-		bdm_assert(alpha.length()==beta.length(), "eBeta:: alpha and beta length do not match");
-		dim = alpha.length();
-	}
-
-	void to_setting ( Setting &set ) const{
-		UI::save(alpha, set, "alpha");
-		UI::save(beta, set, "beta");
-	}
+public:
+    //!sufficient statistics
+    vec alpha;
+    //!sufficient statistics
+    vec beta;
+public:
+    //!\name Constructors
+    //!@{
+
+    eBeta () : eEF () {};
+    eBeta ( const eBeta &B0 ) : eEF (), alpha(B0.alpha),beta(B0.beta) {
+        validate();
+    };
+    //!@}
+
+    //! using sampling procedure from wikipedia
+    vec sample() const {
+        vec y ( beta.length() ); // all vectors
+        for ( int i = 0; i < beta.length(); i++ ) {
+            GamRNG.setup ( alpha ( i ), 1 );
+#pragma omp critical
+            double Ga = GamRNG();
+
+            GamRNG.setup ( beta ( i ), 1 );
+#pragma omp critical
+            double Gb = GamRNG();
+
+            y ( i ) = Ga/(Ga+Gb);
+        }
+        return y;
+    }
+
+    vec mean() const {
+        return elem_div(alpha, alpha + beta); // dot-division
+    };
+    vec variance() const {
+        vec apb=alpha+beta;
+        return elem_div (elem_mult ( alpha, beta) ,
+                         elem_mult ( elem_mult(apb,apb), apb+1 ) );
+    }
+    //! In this instance, val is ...
+    double evallog_nn ( const vec &val ) const {
+        double tmp;
+        tmp = ( alpha - 1 ) * log ( val ) + (beta-1)*log(1-val);
+        return tmp;
+    }
+
+    double lognc () const {
+        double lgb = 0.0;
+        for ( int i = 0; i < beta.length(); i++ ) {
+            lgb += -lgamma ( alpha(i)+beta(i) ) + lgamma(alpha(i)) + lgamma(beta(i));
+        }
+        return lgb;
+    }
+
+    /*! Create object from the following structure
+
+    \code
+    class = 'eBeta';
+    alpha = [...];           % vector parameter alpha
+    beta  = [...];           % vector parameter beta of the same length as alpha
+    \endcode
+
+    Class does not call bdm::eEF::from_setting
+    */
+    void from_setting ( const Setting &set ) {
+        UI::get(alpha, set, "alpha", UI::compulsory);
+        UI::get(beta, set, "beta", UI::compulsory);
+    }
+
+    void validate() {
+        bdm_assert(alpha.length()==beta.length(), "eBeta:: alpha and beta length do not match");
+        dim = alpha.length();
+    }
+
+    void to_setting ( Setting &set ) const {
+        UI::save(alpha, set, "alpha");
+        UI::save(beta, set, "beta");
+    }
 };
 UIREGISTER ( eBeta );
@@ -660,30 +672,30 @@
 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 );
-	void 	to_setting  (Setting  &set) const;
-	void validate();
+    //! 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 );
+    void 	to_setting  (Setting  &set) const;
+    void validate();
 };
 UIREGISTER ( mDirich );
@@ -691,5 +703,5 @@
 /*! \brief Random Walk with vector Beta distribution
 Using simple assignment
-\f{eqnarray*} 
+\f{eqnarray*}
 \alpha & = & rvc / k + \beta_c \\
 \beta & = &(1-rvc) / k + \beta_c \\
@@ -702,39 +714,39 @@
 By default is it set to 0.1;
 */
-class mBeta: public pdf_internal<eBeta>{
-	//! vector of constants \f$ k \f$ of the random walk 
-	vec k;
-	//! stabilizing coefficient \f$ \beta_c \f$
-	vec betac;
-
-	public:
-	void condition ( const vec &val ) {
-		this->iepdf.alpha =  elem_div(val , k) + betac;
-		this->iepdf.beta =  elem_div (1-val , k) + betac;
-	};
-
-	/*! Create Beta random walk
-	\f[ f(rv|rvc) = \prod Beta(rvc,k) \f]
-	from structure
-	\code
-	class = 'mBeta';
-	k = 1;                      // multiplicative constant k
-	--- optional ---
-	rv = RV({'name'},size)      // description of RV
-	beta  = [];                 // initial value of beta
-	betac = [];                 // initial value of beta stabilizing constant
-	\endcode
-	*/
-	void from_setting ( const Setting &set );
-	
-	void to_setting  (Setting  &set) const;
-
-	void validate(){			
-		pdf_internal<eBeta>::validate();
-		bdm_assert(betac.length()==dimension(),"Incomaptible betac");
-		bdm_assert(k.length()==dimension(),"Incomaptible k");
-		dimc = iepdf.dimension();
-	}
-	//! 
+class mBeta: public pdf_internal<eBeta> {
+    //! vector of constants \f$ k \f$ of the random walk
+    vec k;
+    //! stabilizing coefficient \f$ \beta_c \f$
+    vec betac;
+
+public:
+    void condition ( const vec &val ) {
+        this->iepdf.alpha =  elem_div(val , k) + betac;
+        this->iepdf.beta =  elem_div (1-val , k) + betac;
+    };
+
+    /*! Create Beta random walk
+    \f[ f(rv|rvc) = \prod Beta(rvc,k) \f]
+    from structure
+    \code
+    class = 'mBeta';
+    k = 1;                      // multiplicative constant k
+    --- optional ---
+    rv = RV({'name'},size)      // description of RV
+    beta  = [];                 // initial value of beta
+    betac = [];                 // initial value of beta stabilizing constant
+    \endcode
+    */
+    void from_setting ( const Setting &set );
+
+    void to_setting  (Setting  &set) const;
+
+    void validate() {
+        pdf_internal<eBeta>::validate();
+        bdm_assert(betac.length()==dimension(),"Incomaptible betac");
+        bdm_assert(k.length()==dimension(),"Incomaptible k");
+        dimc = iepdf.dimension();
+    }
+    //!
 };
 UIREGISTER(mBeta);
@@ -743,51 +755,51 @@
 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 );
-
-	double logpred ( const vec &yt ) const;
-
-	void flatten ( const BMEF* B , double weight);
-
-	//! return correctly typed posterior (covariant return)
-	const eDirich& posterior() const {
-		return est;
-	};
-	//! constructor function
-	void set_parameters ( const vec &beta0 ) {
-		est.set_parameters ( beta0 );
-		est.validate();
-		if ( evalll ) {
-			last_lognc = est.lognc();
-		}
-	}
-
-	void to_setting ( Setting &set ) const {
-		BMEF::to_setting ( set );
-		UI::save( &est, set, "prior" );
-	}
-	void from_setting (const Setting &set )  {
-		BMEF::from_setting ( set );
-		UI::get( est, set, "prior" );
-	}
+    //! 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 );
+
+    double logpred ( const vec &yt ) const;
+
+    void flatten ( const BMEF* B , double weight);
+
+    //! return correctly typed posterior (covariant return)
+    const eDirich& posterior() const {
+        return est;
+    };
+    //! constructor function
+    void set_parameters ( const vec &beta0 ) {
+        est.set_parameters ( beta0 );
+        est.validate();
+        if ( evalll ) {
+            last_lognc = est.lognc();
+        }
+    }
+
+    void to_setting ( Setting &set ) const {
+        BMEF::to_setting ( set );
+        UI::save( &est, set, "prior" );
+    }
+    void from_setting (const Setting &set )  {
+        BMEF::from_setting ( set );
+        UI::get( est, set, "prior" );
+    }
 };
 UIREGISTER( multiBM );
@@ -804,54 +816,54 @@
 class egamma : public eEF {
 protected:
-	//! Vector \f$\alpha\f$
-	vec alpha;
-	//! Vector \f$\beta\f$
-	vec beta;
+    //! 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 );
-		validate();
-	};
-	void set_parameters ( const vec &a, const vec &b ) {
-		alpha = a, beta = b;
-	};
-	//!@}
-
-	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 object from the following structure
-
-	\code
-	class = 'egamma';
-	alpha = [...];         % vector of alpha
-	beta = [...];          % vector of beta
-	--- inherited fields ---
-	bdm::eEF::from_setting
-	\endcode
-	 fulfilling formula \f[ f(rv|rvc) = \Gamma(\alpha, \beta) \f]
-	*/
-	void from_setting ( const Setting &set );
-
-	void to_setting ( Setting &set ) const;
-	void validate(); 
+    //! \name Constructors
+    //!@{
+    egamma () : eEF (), alpha ( 0 ), beta ( 0 ) {};
+    egamma ( const vec &a, const vec &b ) {
+        set_parameters ( a, b );
+        validate();
+    };
+    void set_parameters ( const vec &a, const vec &b ) {
+        alpha = a, beta = b;
+    };
+    //!@}
+
+    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 object from the following structure
+
+    \code
+    class = 'egamma';
+    alpha = [...];         % vector of alpha
+    beta = [...];          % vector of beta
+    --- inherited fields ---
+    bdm::eEF::from_setting
+    \endcode
+     fulfilling formula \f[ f(rv|rvc) = \Gamma(\alpha, \beta) \f]
+    */
+    void from_setting ( const Setting &set );
+
+    void to_setting ( Setting &set ) const;
+    void validate();
 };
 UIREGISTER ( egamma );
@@ -877,20 +889,20 @@
 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 );
-	}
+    //! \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 );
+    }
 };
 /*
@@ -914,60 +926,60 @@
 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 );
-	}
-	//!@}
-
-	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 );
+    }
+    //!@}
+
+    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 );
-	void 	to_setting  (Setting  &set) const;
-	void validate();
+        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 );
+    void 	to_setting  (Setting  &set) const;
+    void validate();
 };
 UIREGISTER ( euni );
@@ -975,30 +987,30 @@
 //! Uniform density with conditional mean value
 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 );
-	}
-	void 	to_setting  (Setting  &set) const {
-		pdf::to_setting ( set ); 
-		UI::save( iepdf.mean(), set, "delta");
-		UI::save(mean, set, "mean");
-	}
-	void validate(){
-		pdf_internal<euni>::validate();
-		dimc = mean->dimensionc();
-		
-	}
+    //! 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 );
+    }
+    void 	to_setting  (Setting  &set) const {
+        pdf::to_setting ( set );
+        UI::save( iepdf.mean(), set, "delta");
+        UI::save(mean, set, "mean");
+    }
+    void validate() {
+        pdf_internal<euni>::validate();
+        dimc = mean->dimensionc();
+
+    }
 
 };
@@ -1012,91 +1024,91 @@
 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;
+    //! 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 );
-		validate();
-	}
-
-	//! 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;
-	}
-
-	//!@}
-	//! 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;
+    //! \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 );
+        validate();
+    }
+
+    //! 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;
+    }
+
+    //!@}
+    //! 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
-	R     = [];                  // square matrix of appropriate dimension
-	--- optional ---
-	const = zeros(A.rows);       // vector of constant term
-	\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::optional);
-		mat R0;
-		UI::get ( R0, set, "R", UI::compulsory );
-		set_parameters ( A, mu_const, R0 );
-	}
-
-void to_setting (Setting &set) const {
-		pdf::to_setting(set);
-		UI::save ( A, set, "A");
-		UI::save ( mu_const, set, "const");
-		UI::save ( _R(), set, "R");
-	}
-
-void validate() {
-		pdf_internal<TEpdf<sq_T> >::validate();
-		if (mu_const.length()==0) { // default in from_setting
-			mu_const=zeros(A.rows());
-		}
-		bdm_assert ( A.rows() == mu_const.length(), "mlnorm: A vs. mu mismatch" );
-		bdm_assert ( A.rows() == _R().rows(), "mlnorm: A vs. R mismatch" );
-		this->dimc = A.cols();
-
-	}
+    }
+
+    //!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
+    R     = [];                  // square matrix of appropriate dimension
+    --- optional ---
+    const = zeros(A.rows);       // vector of constant term
+    \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::optional);
+        mat R0;
+        UI::get ( R0, set, "R", UI::compulsory );
+        set_parameters ( A, mu_const, R0 );
+    }
+
+    void to_setting (Setting &set) const {
+        pdf::to_setting(set);
+        UI::save ( A, set, "A");
+        UI::save ( mu_const, set, "const");
+        UI::save ( _R(), set, "R");
+    }
+
+    void validate() {
+        pdf_internal<TEpdf<sq_T> >::validate();
+        if (mu_const.length()==0) { // default in from_setting
+            mu_const=zeros(A.rows());
+        }
+        bdm_assert ( A.rows() == mu_const.length(), "mlnorm: A vs. mu mismatch" );
+        bdm_assert ( A.rows() == _R().rows(), "mlnorm: A vs. R mismatch" );
+        this->dimc = A.cols();
+
+    }
 };
 UIREGISTER2 ( mlnorm, ldmat );
@@ -1107,5 +1119,5 @@
 SHAREDPTR2 ( mlnorm, chmat );
 
-//! \class mlgauss 
+//! \class mlgauss
 //!\brief Normal distribution with linear function of mean value. Same as mlnorm<fsqmat>.
 typedef mlnorm<fsqmat> mlgauss;
@@ -1117,64 +1129,64 @@
 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 to_setting  (Setting  &set) const {
-		UI::save( g,set, "g");
-		UI::save(this->iepdf._R().to_mat(),set, "R");
-	
-	}
-
-
-
-void validate() {
-		this->iepdf.validate();
-		bdm_assert ( g->dimension() == this->iepdf.dimension(), "incompatible function" );
-		this->dim = g->dimension();
-		this->dimc = g->dimensionc();
-		this->iepdf.validate();
-	}
+    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 to_setting  (Setting  &set) const {
+        UI::save( g,set, "g");
+        UI::save(this->iepdf._R().to_mat(),set, "R");
+
+    }
+
+
+
+    void validate() {
+        this->iepdf.validate();
+        bdm_assert ( g->dimension() == this->iepdf.dimension(), "incompatible function" );
+        this->dim = g->dimension();
+        this->dimc = g->dimensionc();
+        this->iepdf.validate();
+    }
 
 };
@@ -1194,30 +1206,30 @@
 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 ); 
-
-	void validate() {
-		mlnorm<ldmat, enorm>::validate();
-		bdm_assert ( A.rows() == mu_const.length(), "mlstudent: A vs. mu mismatch" );
-		bdm_assert ( _R.rows() == A.rows(), "mlstudent: A vs. R mismatch" );
-
-	}
+    //! 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 );
+
+    void validate() {
+        mlnorm<ldmat, enorm>::validate();
+        bdm_assert ( A.rows() == mu_const.length(), "mlstudent: A vs. mu mismatch" );
+        bdm_assert ( _R.rows() == A.rows(), "mlstudent: A vs. R mismatch" );
+
+    }
 };
 
@@ -1234,36 +1246,36 @@
 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 );
-	void 	to_setting  (Setting  &set) const;
-	void validate();
+    //! 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 );
+    void 	to_setting  (Setting  &set) const;
+    void validate();
 };
 UIREGISTER ( mgamma );
@@ -1281,43 +1293,43 @@
 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*/ );		
-	};
-
-	void 	validate (){
-		pdf_internal<eigamma>::validate();
-		dimc = dimension();
-};
-
-	void condition ( const vec &val ) {
-		_beta = elem_mult ( val, ( _alpha - 1.0 ) );
-	};
+    //! 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*/ );
+    };
+
+    void 	validate () {
+        pdf_internal<eigamma>::validate();
+        dimc = dimension();
+    };
+
+    void condition ( const vec &val ) {
+        _beta = elem_mult ( val, ( _alpha - 1.0 ) );
+    };
 };
 
@@ -1336,27 +1348,27 @@
 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;	
-	};
-
-	void 	validate (){
-		mgamma::validate();
-		dimc = dimension();
-	};
-
-	void condition ( const vec &val ) {
-		vec mean = elem_mult ( refl, pow ( val, l ) );
-		_beta = k / mean;
-	};
+    //! 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;
+    };
+
+    void 	validate () {
+        mgamma::validate();
+        dimc = dimension();
+    };
+
+    void condition ( const vec &val ) {
+        vec mean = elem_mult ( refl, pow ( val, l ) );
+        _beta = k / mean;
+    };
 };
 
@@ -1376,45 +1388,45 @@
 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;
-	};
-
-	void validate(){
-		migamma::validate();
-		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 );
-
-	void to_setting  (Setting  &set) const; 
+    //! 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;
+    };
+
+    void validate() {
+        migamma::validate();
+        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 );
+
+    void to_setting  (Setting  &set) const;
 };
 
@@ -1435,11 +1447,11 @@
 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 );
-	};
+    vec sample() const {
+        return exp ( enorm<ldmat>::sample() );
+    };
+    vec mean() const {
+        vec var = enorm<ldmat>::variance();
+        return exp ( mu - 0.5*var );
+    };
 
 };
@@ -1453,45 +1465,45 @@
 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 ) );
-	};
-	
-	void validate(){
-		pdf_internal<elognorm>::validate();
-		dimc = iepdf.dimension();
-	}
-
-	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 );
-	
-	void to_setting  (Setting  &set) const; 
+    //! 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 ) );
+    };
+
+    void validate() {
+        pdf_internal<elognorm>::validate();
+        dimc = iepdf.dimension();
+    }
+
+    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 );
+
+    void to_setting  (Setting  &set) const;
 };
 
@@ -1504,74 +1516,74 @@
 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();	
-	}
-	//! Set internal structures
-	void set_parameters ( const chmat &Y0, const double delta0 ) {
-		Y = Y0;
-		delta = delta0;
-		p = Y.rows();
-		}
-	
-	virtual void validate (){
-		epdf::validate();
-		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
+    //! 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();
+    }
+    //! Set internal structures
+    void set_parameters ( const chmat &Y0, const double delta0 ) {
+        Y = Y0;
+        delta = delta0;
+        p = Y.rows();
+    }
+
+    virtual void validate () {
+        epdf::validate();
+        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() );
-		}
-		//do the rest
-		for ( int i = 0; i < p; i++ ) {
-			for ( int j = i + 1; j < p; j++ ) {
+            X ( i, i ) = sqrt ( GamRNG() );
+        }
+        //do the rest
+        for ( int i = 0; i < p; i++ ) {
+            for ( int j = i + 1; j < p; j++ ) {
 #pragma omp critical
-				X ( i, j ) = NorRNG.sample();
-			}
-		}
-		return X*Y._Ch();// return upper triangular part of the decomposition
-	}
-
-	vec sample () const {
-		return vec ( sample_mat()._data(), p*p );
-	}
-
-	virtual vec mean() const NOT_IMPLEMENTED(0);
-
-	//! return expected variance (not covariance!)
-	virtual vec variance() const NOT_IMPLEMENTED(0);
-
-	virtual double evallog ( const vec &val ) const NOT_IMPLEMENTED(0);
-
-	//! 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;
-	}
+                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 );
+    }
+
+    virtual vec mean() const NOT_IMPLEMENTED(0);
+
+    //! return expected variance (not covariance!)
+    virtual vec variance() const NOT_IMPLEMENTED(0);
+
+    virtual double evallog ( const vec &val ) const NOT_IMPLEMENTED(0);
+
+    //! 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;
+    }
 };
 
@@ -1581,68 +1593,68 @@
 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 );
-		p = Y0.rows();
-	}
-	
-	virtual void 	validate (){
-		epdf::validate();
-		W.validate();
-		dim = W.dimension();
-	}
-	
-	
-	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
+    //! 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 );
+        p = Y0.rows();
+    }
+
+    virtual void 	validate () {
+        epdf::validate();
+        W.validate();
+        dim = W.dimension();
+    }
+
+
+    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;
-	};
-
-	virtual vec mean() const NOT_IMPLEMENTED(0);
-
-	//! return expected variance (not covariance!)
-	virtual vec variance() const NOT_IMPLEMENTED(0);
+        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;
+    };
+
+    virtual vec mean() const NOT_IMPLEMENTED(0);
+
+    //! return expected variance (not covariance!)
+    virtual vec variance() const NOT_IMPLEMENTED(0);
 };
 
@@ -1650,40 +1662,40 @@
 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 );
-	};
-	
-	void validate(){
-		pdf_internal<eiWishartCh>::validate();
-		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 );
-	}
+    //!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 );
+    };
+
+    void validate() {
+        pdf_internal<eiWishartCh>::validate();
+        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 );
+    }
 };
 
@@ -1691,5 +1703,5 @@
 enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
 
-//! Shortcut for multinomial.sample(int n). Various simplifications are allowed see RESAMPLING_METHOD 
+//! Shortcut for multinomial.sample(int n). Various simplifications are allowed see RESAMPLING_METHOD
 void resample(const vec &w, ivec &ind, RESAMPLING_METHOD=SYSTEMATIC);
 
@@ -1700,98 +1712,98 @@
 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 ) {
-		n = Av.size();
-		w = 1 / n * ones ( n );
-		samples = Av;
-	};
-	virtual void 	validate ();
-	//! 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.
-	void resample ( RESAMPLING_METHOD method = SYSTEMATIC );
-
-	//! inherited operation : NOT implemented
-	vec sample() const NOT_IMPLEMENTED(0);
-
-	//! inherited operation : NOT implemented
-	double evallog ( const vec &val ) const NOT_IMPLEMENTED(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;
-
-	void to_setting ( Setting &set ) const;
-	
-	/*! Create object from the following structure
-
-	\code
-	class = 'eEmp';
-	samples = [...];               % array of samples 
-	w = [...];					   % weights of samples stored in vector
-	--- inherited fields ---
-	bdm::epdf::from_setting	
-	\endcode	
-	*/
-	void from_setting ( const Setting &set );
+    //! 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 ) {
+        n = Av.size();
+        w = 1 / n * ones ( n );
+        samples = Av;
+    };
+    virtual void 	validate ();
+    //! 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.
+    void resample ( RESAMPLING_METHOD method = SYSTEMATIC );
+
+    //! inherited operation : NOT implemented
+    vec sample() const NOT_IMPLEMENTED(0);
+
+    //! inherited operation : NOT implemented
+    double evallog ( const vec &val ) const NOT_IMPLEMENTED(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;
+
+    void to_setting ( Setting &set ) const;
+
+    /*! Create object from the following structure
+
+    \code
+    class = 'eEmp';
+    samples = [...];               % array of samples
+    w = [...];					   % weights of samples stored in vector
+    --- inherited fields ---
+    bdm::epdf::from_setting
+    \endcode
+    */
+    void from_setting ( const Setting &set );
 };
 UIREGISTER(eEmp);
@@ -1803,31 +1815,31 @@
 void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) {
 //Fixme test dimensions of mu0 and R0;
-	mu = mu0;
-	R = R0;
-	validate();
+    mu = mu0;
+    R = R0;
+    validate();
 };
 
 template<class sq_T>
 void enorm<sq_T>::from_setting ( const Setting &set ) {
-	epdf::from_setting ( set ); //reads rv
-
-	UI::get ( mu, set, "mu", UI::compulsory );
-	mat Rtmp;// necessary for conversion
-	UI::get ( Rtmp, set, "R", UI::compulsory );
-	R = Rtmp; // conversion
+    epdf::from_setting ( set ); //reads rv
+
+    UI::get ( mu, set, "mu", UI::compulsory );
+    mat Rtmp;// necessary for conversion
+    UI::get ( Rtmp, set, "R", UI::compulsory );
+    R = Rtmp; // conversion
 }
 
 template<class sq_T>
 void enorm<sq_T>::validate() {
-		eEF::validate();
-		bdm_assert ( mu.length() == R.rows(), "mu and R parameters do not match" );
-		dim = mu.length();
-	}
+    eEF::validate();
+    bdm_assert ( mu.length() == R.rows(), "mu and R parameters do not match" );
+    dim = mu.length();
+}
 
 template<class sq_T>
 void enorm<sq_T>::to_setting ( Setting &set ) const {
-	epdf::to_setting ( set ); //reads rv	
-	UI::save ( mu, set, "mu");
-	UI::save ( R.to_mat(), set, "R");
+    epdf::to_setting ( set ); //reads rv
+    UI::save ( mu, set, "mu");
+    UI::save ( R.to_mat(), set, "R");
 }
 
@@ -1836,5 +1848,5 @@
 template<class sq_T>
 void enorm<sq_T>::dupdate ( mat &v, double nu ) {
-	//
+    //
 };
 
@@ -1846,11 +1858,11 @@
 template<class sq_T>
 vec enorm<sq_T>::sample() const {
-	vec x ( dim );
+    vec x ( dim );
 #pragma omp critical
-	NorRNG.sample_vector ( dim, x );
-	vec smp = R.sqrt_mult ( x );
-
-	smp += mu;
-	return smp;
+    NorRNG.sample_vector ( dim, x );
+    vec smp = R.sqrt_mult ( x );
+
+    smp += mu;
+    return smp;
 };
 
@@ -1865,14 +1877,14 @@
 template<class sq_T>
 double enorm<sq_T>::evallog_nn ( const vec &val ) const {
-	// 1.83787706640935 = log(2pi)
-	double tmp = -0.5 * ( R.invqform ( mu - val ) );// - lognc();
-	return  tmp;
+    // 1.83787706640935 = log(2pi)
+    double tmp = -0.5 * ( R.invqform ( mu - val ) );// - lognc();
+    return  tmp;
 };
 
 template<class sq_T>
 inline double enorm<sq_T>::lognc () const {
-	// 1.83787706640935 = log(2pi)
-	double tmp = 0.5 * ( R.cols() * 1.83787706640935 + R.logdet() );
-	return tmp;
+    // 1.83787706640935 = log(2pi)
+    double tmp = 0.5 * ( R.cols() * 1.83787706640935 + R.logdet() );
+    return tmp;
 };
 
@@ -1906,78 +1918,94 @@
 template<class sq_T>
 shared_ptr<epdf> enorm<sq_T>::marginal ( const RV &rvn ) const {
-	enorm<sq_T> *tmp = new enorm<sq_T> ();
-	shared_ptr<epdf> narrow ( tmp );
-	marginal ( rvn, *tmp );
-	return narrow;
+    enorm<sq_T> *tmp = new enorm<sq_T> ();
+    shared_ptr<epdf> narrow ( tmp );
+    marginal ( rvn, *tmp );
+    return narrow;
 }
 
 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
-
-	target.set_rv ( rvn );
-	target.set_parameters ( mu ( irvn ), Rn );
+    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 );
 }
 
 template<class sq_T>
 shared_ptr<pdf> enorm<sq_T>::condition ( const RV &rvn ) const {
-	mlnorm<sq_T> *tmp = new mlnorm<sq_T> ();
-	shared_ptr<pdf> narrow ( tmp );
-	condition ( rvn, *tmp );
-	return narrow;
+    mlnorm<sq_T> *tmp = new mlnorm<sq_T> ();
+    shared_ptr<pdf> narrow ( tmp );
+    condition ( rvn, *tmp );
+    return narrow;
 }
 
 template<class sq_T>
 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" );
-	//Permutation vector of the new R
-	ivec irvn = rvn.dataind ( rv );
-	ivec irvc = rvc.dataind ( rv );
-	ivec perm = concat ( irvn , irvc );
-	sq_T Rn ( R, perm );
-
-	//fixme - could this be done in general for all sq_T?
-	mat S = Rn.to_mat();
-	//fixme
-	int n = rvn._dsize() - 1;
-	int end = R.rows() - 1;
-	mat S11 = S.get ( 0, n, 0, n );
-	mat S12 = S.get ( 0, n , rvn._dsize(), end );
-	mat S22 = S.get ( rvn._dsize(), end, rvn._dsize(), end );
-
-	vec mu1 = mu ( irvn );
-	vec mu2 = mu ( irvc );
-	mat A = S12 * inv ( S22 );
-	sq_T R_n ( S11 - A *S12.T() );
-
-	uptarget.set_rv ( rvn );
-	uptarget.set_rvc ( rvc );
-	uptarget.set_parameters ( A, mu1 - A*mu2, R_n );
-	uptarget.validate();
+    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" );
+    //Permutation vector of the new R
+    ivec irvn = rvn.dataind ( rv );
+    ivec irvc = rvc.dataind ( rv );
+    ivec perm = concat ( irvn , irvc );
+    sq_T Rn ( R, perm );
+
+    //fixme - could this be done in general for all sq_T?
+    mat S = Rn.to_mat();
+    //fixme
+    int n = rvn._dsize() - 1;
+    int end = R.rows() - 1;
+    mat S11 = S.get ( 0, n, 0, n );
+    mat S12 = S.get ( 0, n , rvn._dsize(), end );
+    mat S22 = S.get ( rvn._dsize(), end, rvn._dsize(), end );
+
+    vec mu1 = mu ( irvn );
+    vec mu2 = mu ( irvc );
+    mat A = S12 * inv ( S22 );
+    sq_T R_n ( S11 - A *S12.T() );
+
+    uptarget.set_rv ( rvn );
+    uptarget.set_rvc ( rvc );
+    uptarget.set_parameters ( A, mu1 - A*mu2, R_n );
+    uptarget.validate();
 }
 
 /*! \brief Dirac delta function distribution */
-class dirac: public epdf{
-	public: 
-		vec point;
-	public:
-		double evallog (const vec &dt) const {return -inf;}
-		vec mean () const {return point;}
-		vec variance () const {return zeros(point.length());}
-		void qbounds ( vec &lb, vec &ub, double percentage = 0.95 ) const { lb = point; ub = point;}
-		//! access
-		const vec& _point() {return point;}
-		void set_point(const vec& p){point =p; dim=p.length();}
-		vec sample() const {return point;}
-	};
+class dirac: public epdf {
+public:
+    vec point;
+public:
+    double evallog (const vec &dt) const {
+        return -inf;
+    }
+    vec mean () const {
+        return point;
+    }
+    vec variance () const {
+        return zeros(point.length());
+    }
+    void qbounds ( vec &lb, vec &ub, double percentage = 0.95 ) const {
+        lb = point;
+        ub = point;
+    }
+    //! access
+    const vec& _point() {
+        return point;
+    }
+    void set_point(const vec& p) {
+        point =p;
+        dim=p.length();
+    }
+    vec sample() const {
+        return point;
+    }
+};
 
 
@@ -1986,11 +2014,11 @@
 template<class sq_T>
 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 );
+    g = g0;
+    this->iepdf.set_parameters ( zeros ( g->dimension() ), R0 );
 }
 
 template<class sq_T>
 void mgnorm<sq_T >::condition ( const vec &cond ) {
-	this->iepdf._mu() = g->eval ( cond );
+    this->iepdf._mu() = g->eval ( cond );
 };
 
@@ -1998,8 +2026,8 @@
 template<class sq_T>
 std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_T> &ml ) {
-	os << "A:" << ml.A << endl;
-	os << "mu:" << ml.mu_const << endl;
-	os << "R:" << ml._R() << endl;
-	return os;
+    os << "A:" << ml.A << endl;
+    os << "mu:" << ml.mu_const << endl;
+    os << "R:" << ml._R() << endl;
+    return os;
 };
 
