/*!
  \file
  \brief Probability distributions for Exponential Family models.
  \author Vaclav Smidl.

  -----------------------------------
  BDM++ - C++ library for Bayesian Decision Making under Uncertainty

  Using IT++ for numerical operations
  -----------------------------------
*/

#ifndef EF_H
#define EF_H


#include "../shared_ptr.h"
#include "../base/bdmbase.h"
#include "../math/chmat.h"

namespace bdm {


//! Global Uniform_RNG
extern Uniform_RNG UniRNG;
//! Global Normal_RNG
extern Normal_RNG NorRNG;
//! Global Gamma_RNG
extern Gamma_RNG GamRNG;

/*!
* \brief General conjugate exponential family posterior density.

* More?...
*/

class eEF : public epdf {
public:
//	eEF() :epdf() {};
	//! default constructor
	eEF () : epdf () {};
	//! logarithm of the normalizing constant, \f$\mathcal{I}\f$
	virtual double lognc() const = 0;

	//!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;
};


//! 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 ) 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 ) {};
	//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 ) NOT_IMPLEMENTED_VOID;


	void to_setting ( Setting &set ) const
	{			
		BM::to_setting( set );
		UI::save(frg, set, "frg");
	} 

	void from_setting( const Setting &set) {
		BM::from_setting(set);
		if ( !UI::get ( frg, set, "frg" ) )
			frg = 1.0;
	}

	void validate() {
		BM::validate();	
	}

};

/*! Dirac delta density with predefined transformation

Density of the type:\f[ f(x_t | y_t) = \delta (x_t - g(y_t)) \f]
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){
			pdf::from_setting(set);
			g=UI::build<fnc>(set,"g",UI::compulsory);
			validate();
		}
		void to_setting(Setting &set) const{
			pdf::to_setting(set);
			UI::save(g.get(), set, "g");
		}
		void validate() {
			dim = g->dimension();
			dimc = g->dimensionc();
		}
};
UIREGISTER(mgdirac);


template<class sq_T, template <typename> class TEpdf>
class mlnorm;

/*!
* \brief Gaussian density with positive definite (decomposed) covariance matrix.

* More?...
*/
template<class sq_T>
class enorm : public eEF {
protected:
	//! mean value
	vec mu;
	//! Covariance matrix in decomposed form
	sq_T R;
public:
	//!\name Constructors
	//!@{

	enorm () : eEF (), mu (), R () {};
	enorm ( const vec &mu, const sq_T &R ) {
		set_parameters ( mu, R );
	}
	void set_parameters ( const vec &mu, const sq_T &R );
	/*! 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() {
		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 );

	//! 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;
	}
	//!@}

};
UIREGISTER2 ( enorm, chmat );
SHAREDPTR2 ( enorm, chmat );
UIREGISTER2 ( enorm, ldmat );
SHAREDPTR2 ( enorm, ldmat );
UIREGISTER2 ( enorm, fsqmat );
SHAREDPTR2 ( enorm, fsqmat );

typedef enorm<ldmat> egauss;
UIREGISTER(egauss);


//forward declaration
class mstudent;

/*! distribution of multivariate Student t density

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() {
				dim = H.rows();
			}
};
UIREGISTER2(estudent,fsqmat);
UIREGISTER2(estudent,ldmat);
UIREGISTER2(estudent,chmat);

/*!
* \brief Gauss-inverse-Wishart density stored in LD form

* For \f$p\f$-variate densities, given rv.count() should be \f$p\times\f$ V.rows().
*
*/
class egiw : public eEF {
	//! \var log_level_enums logvariance
	//! TODO DOPLNIT

	//! \var log_level_enums logmean
	//! TODO DOPLNIT
	LOG_LEVEL(egiw,logmean, logvariance);

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 );
		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;
	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     = [];               // 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
	log_level = 'tri';		  // set the level of logged details
	\endcode

	\sa log_level_enums
	*/
	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 );
			validate();

		} else {
			set_parameters ( dimx, V, nu );
			validate();
		}
	}

	void to_setting ( Setting& set ) const {
		epdf::to_setting ( set );
		UI::save ( dimx, set, "dimx" );
		UI::save ( V.to_mat(), set, "V" );
		UI::save ( nu, set, "nu" );
	};

	void validate() {
	dim = dimx * ( dimx + nPsi ); 	
	
	    // check sizes, rvs etc.
		// also check if RV are meaningful!!!
		// meaningful =  rv for theta  and rv for r are split!
	}
	void log_register ( bdm::logger& L, const string& prefix );

	void log_write() const;
	//!@}
};
UIREGISTER ( egiw );
SHAREDPTR ( egiw );

/*! \brief Dirichlet posterior density

Continuous Dirichlet density of \f$n\f$-dimensional variable \f$x\f$
\f[
f(x|\beta) = \frac{\Gamma[\gamma]}{\prod_{i=1}^{n}\Gamma(\beta_i)} \prod_{i=1}^{n}x_i^{\beta_i-1}
\f]
where \f$\gamma=\sum_i \beta_i\f$.
*/
class eDirich: public eEF {
protected:
	//!sufficient statistics
	vec beta;
public:
	//!\name Constructors
	//!@{

	eDirich () : eEF () {};
	eDirich ( const eDirich &D0 ) : eEF () {
		set_parameters ( D0.beta );
		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;
	}
	/*! 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();
	}

	void to_setting ( Setting &set ) const
	{			
		eEF::to_setting( set );
		UI::save( beta, set, "beta" );
	} 
};
UIREGISTER ( eDirich );

/*! Random Walk on Dirichlet
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;
*/

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 );

//! \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 );

	double logpred ( const vec &yt ) const;

	void flatten ( const BMEF* B );

	//! 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" );
	}
};
UIREGISTER( multiBM );

/*!
 \brief Gamma posterior density

 Multivariate Gamma density as product of independent univariate densities.
 \f[
 f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i)
 \f]
*/

class egamma : public eEF {
protected:
	//! Vector \f$\alpha\f$
	vec alpha;
	//! Vector \f$\beta\f$
	vec beta;
public :
	//! \name Constructors
	//!@{
	egamma () : eEF (), alpha ( 0 ), beta ( 0 ) {};
	egamma ( const vec &a, const vec &b ) {
		set_parameters ( a, b );
		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 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 );

/*!
 \brief Inverse-Gamma posterior density

 Multivariate inverse-Gamma density as product of independent univariate densities.
 \f[
 f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i)
 \f]

Vector \f$\beta\f$ has different meaning (in fact it is 1/beta as used in definition of iG)

 Inverse Gamma can be converted to Gamma using \f[
 x\sim iG(a,b) => 1/x\sim G(a,1/b)
\f]
This relation is used in sampling.
 */

class eigamma : public egamma {
protected:
public :
	//! \name Constructors
	//! All constructors are inherited
	//!@{
	//!@}

	vec sample() const {
		return 1.0 / egamma::sample();
	};
	//! Returns poiter to alpha and beta. Potentially dangerous: use with care!
	vec mean() const {
		return elem_div ( beta, alpha - 1 );
	}
	vec variance() const {
		vec mea = mean();
		return elem_div ( elem_mult ( mea, mea ), alpha - 2 );
	}
};
/*
//! Weighted mixture of epdfs with external owned components.
class emix : public epdf {
protected:
	int n;
	vec &w;
	Array<epdf*> Coms;
public:
//! Default constructor
	emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
	void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
	vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
};
*/

//!  Uniform distributed density on a rectangular support

class euni: public epdf {
protected:
//! lower bound on support
	vec low;
//! upper bound on support
	vec high;
//! internal
	vec distance;
//! normalizing coefficients
	double nk;
//! cache of log( \c nk )
	double lnk;
public:
	//! \name Constructors
	//!@{
	euni () : epdf () {}
	euni ( const vec &low0, const vec &high0 ) {
		set_parameters ( low0, high0 );
	}
	void set_parameters ( const vec &low0, const vec &high0 ) {
		distance = high0 - low0;
		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 ) {
		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> {
	//! 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 validate(){
		dimc = mean->dimensionc();
		iepdf.validate();
	}

};
UIREGISTER ( mguni );
/*!
 \brief Normal distributed linear function with linear function of mean value;

 Mean value \f$ \mu=A*\mbox{rvc}+\mu_0 \f$.
*/
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;
//			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;
//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 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();
		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 );
SHAREDPTR2 ( mlnorm, ldmat );
UIREGISTER2 ( mlnorm, fsqmat );
SHAREDPTR2 ( mlnorm, fsqmat );
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:
//			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" );
		this->dim = g->dimension();
		this->dimc = g->dimensionc();
		this->iepdf.validate();
	}
};

UIREGISTER2 ( mgnorm, chmat );
SHAREDPTR2 ( mgnorm, chmat );


/*! (Approximate) Student t density with linear function of mean value

The internal epdf of this class is of the type of a Gaussian (enorm).
However, each conditioning is trying to assure the best possible approximation by taking into account the zeta function. See [] for reference.

Perhaps a moment-matching technique?
*/
class mlstudent : public mlnorm<ldmat, enorm> {
protected:
	//! 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() {
		bdm_assert ( A.rows() == mu_const.length(), "mlstudent: A vs. mu mismatch" );
		bdm_assert ( _R.rows() == A.rows(), "mlstudent: A vs. R mismatch" );

	}
};

/*!
\brief  Gamma random walk

Mean value, \f$\mu\f$, of this density is given by \c rvc .
Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.

The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
*/
class mgamma : public 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 );

/*!
\brief  Inverse-Gamma random walk

Mean value, \f$ \mu \f$, of this density is given by \c rvc .
Standard deviation of the random walk is proportional to one \f$ k \f$-th the mean.
This is achieved by setting \f$ \alpha=\mu/k^2+2 \f$ and \f$ \beta=\mu(\alpha-1)\f$.

The standard deviation of the walk is then: \f$ \mu/\sqrt(k)\f$.
 */
class migamma : public 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 (){
		iepdf.validate();
		dimc = dimension();
};

	void condition ( const vec &val ) {
		_beta = elem_mult ( val, ( _alpha - 1.0 ) );
	};
};


/*!
\brief  Gamma random walk around a fixed point

Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]

Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.

The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
*/
class mgamma_fix : public mgamma {
protected:
	//! parameter l
	double l;
	//! reference vector
	vec refl;
public:
	//! Constructor
	mgamma_fix () : mgamma (), refl () {};
	//! Set value of \c k
	void set_parameters ( double k0 , vec ref0, double l0 ) {
		mgamma::set_parameters ( k0, ref0 );
		refl = pow ( ref0, 1.0 - l0 );
		l = l0;	
	};

	void 	validate (){
		mgamma::validate();
		dimc = dimension();
	};

	void condition ( const vec &val ) {
		vec mean = elem_mult ( refl, pow ( val, l ) );
		_beta = k / mean;
	};
};


/*!
\brief  Inverse-Gamma random walk around a fixed point

Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]

==== Check == vv =
Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.

The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
 */
class migamma_ref : public migamma {
protected:
	//! parameter l
	double l;
	//! reference vector
	vec refl;
public:
	//! Constructor
	migamma_ref () : migamma (), refl () {};
	
	//! Set value of \c k
	void set_parameters ( double k0 , vec ref0, double l0 ) {
		migamma::set_parameters ( ref0.length(), k0 );
		refl = pow ( ref0, 1.0 - l0 );
		l = l0;
	};

	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 );

	// TODO dodelat void to_setting( Setting &set ) const;
};


UIREGISTER ( migamma_ref );
SHAREDPTR ( migamma_ref );

/*! Log-Normal probability density
 only allow diagonal covariances!

Density of the form \f$ \log(x)\sim \mathcal{N}(\mu,\sigma^2) \f$ , i.e.
\f[
x \sim \frac{1}{x\sigma\sqrt{2\pi}}\exp{-\frac{1}{2\sigma^2}(\log(x)-\mu)}
\f]

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 );
	};

};

/*!
\brief  Log-Normal random walk

Mean value, \f$\mu\f$, is...

 */
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(){
		iepdf.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 );

	// 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();	
	}
	//! Set internal structures
	void set_parameters ( const chmat &Y0, const double delta0 ) {
		Y = Y0;
		delta = delta0;
		p = Y.rows();
		}
	
	virtual void 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++ ) {
#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;
	}
};

//! Inverse Wishart on Choleski decomposition
/*! 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 );
		p = Y0.rows();
	}
	
	virtual void 	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);
};

//! 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 );
	};
	
	void validate(){
		iepdf.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 );
	}
};

//! Switch between various resampling methods.
enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };

//! Shortcut for multinomial.sample(int n). Various simplifications are allowed see RESAMPLING_METHOD 
void resample(const vec &w, ivec &ind, RESAMPLING_METHOD=SYSTEMATIC);

/*!
\brief Weighted empirical density

Used e.g. in particle filters.
*/
class eEmp: public epdf {
protected :
	//! Number of particles
	int n;
	//! Sample weights \f$w\f$
	vec w;
	//! Samples \f$x^{(i)}, i=1..n\f$
	Array<vec> samples;
public:
	//! \name Constructors
	//!@{
	eEmp () : epdf (), w (), samples () {};
	//! copy constructor
	eEmp ( const eEmp &e ) : epdf ( e ), w ( e.w ), samples ( e.samples ) {};
	//!@}

	//! Set samples and weights
	void set_statistics ( const vec &w0, const epdf &pdf0 );
	//! Set samples and weights
	void set_statistics ( const epdf &pdf0 , int n ) {
		set_statistics ( ones ( n ) / n, pdf0 );
	};
	//! Set sample
	void set_samples ( const epdf* pdf0 );
	//! Set sample
	void set_parameters ( int n0, bool copy = true ) {
		n = n0;
		w.set_size ( n0, copy );
		samples.set_size ( n0, copy );
	};
	//! Set samples
	void set_parameters ( const Array<vec> &Av ) {
		n = Av.size();
		w = 1 / n * ones ( n );
		samples = Av;
	};
	virtual void 	validate (){
	  bdm_assert (samples.length()==w.length(),"samples and weigths are of different lengths");
	  n = w.length();
	  if (n>0)
		pdf::dim = samples ( 0 ).length();
	}
	//! 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 {
		epdf::to_setting( set );
		UI::save ( samples, set, "samples" );
		UI::save ( w, set, "w" );
	}

	void from_setting ( const Setting &set ) {
		epdf::from_setting( set );
		
		UI::get( samples, set, "samples", UI::compulsory );
		UI::get ( w, set, "w", UI::compulsory );
		validate();
	}

};
UIREGISTER(eEmp);


////////////////////////

template<class sq_T>
void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 ) {
//Fixme test dimensions of mu0 and R0;
	mu = mu0;
	R = R0;
	validate();
};

template<class sq_T>
void enorm<sq_T>::from_setting ( const Setting &set ) {
	epdf::from_setting ( set ); //reads rv

	UI::get ( mu, set, "mu", UI::compulsory );
	mat Rtmp;// necessary for conversion
	UI::get ( Rtmp, set, "R", UI::compulsory );
	R = Rtmp; // conversion
	validate();
}

template<class sq_T>
void enorm<sq_T>::to_setting ( Setting &set ) const {
	epdf::to_setting ( set ); //reads rv	
	UI::save ( mu, set, "mu");
	UI::save ( R.to_mat(), set, "R");
}

template<class sq_T>
void enorm<sq_T>::dupdate ( mat &v, double nu ) {
	//
};

// template<class sq_T>
// void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
// 	//
// };

template<class sq_T>
vec enorm<sq_T>::sample() const {
	vec x ( dim );
#pragma omp critical
	NorRNG.sample_vector ( dim, x );
	vec smp = R.sqrt_mult ( x );

	smp += mu;
	return smp;
};

// template<class sq_T>
// double enorm<sq_T>::eval ( const vec &val ) const {
// 	double pdfl,e;
// 	pdfl = evallog ( val );
// 	e = exp ( pdfl );
// 	return e;
// };

template<class sq_T>
double enorm<sq_T>::evallog_nn ( const vec &val ) const {
	// 1.83787706640935 = log(2pi)
	double tmp = -0.5 * ( R.invqform ( mu - val ) );// - lognc();
	return  tmp;
};

template<class sq_T>
inline double enorm<sq_T>::lognc () const {
	// 1.83787706640935 = log(2pi)
	double tmp = 0.5 * ( R.cols() * 1.83787706640935 + R.logdet() );
	return tmp;
};


// template<class sq_T>
// vec mlnorm<sq_T>::samplecond (const  vec &cond, double &lik ) {
// 	this->condition ( cond );
// 	vec smp = epdf.sample();
// 	lik = epdf.eval ( smp );
// 	return smp;
// }

// template<class sq_T>
// mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) {
// 	int i;
// 	int dim = rv.count();
// 	mat Smp ( dim,n );
// 	vec smp ( dim );
// 	this->condition ( cond );
//
// 	for ( i=0; i<n; i++ ) {
// 		smp = epdf.sample();
// 		lik ( i ) = epdf.eval ( smp );
// 		Smp.set_col ( i ,smp );
// 	}
//
// 	return Smp;
// }


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;
}

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 );
}

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;
}

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();
}

/*! Dirac delta function distribution */
class dirac: public epdf{
	protected: 
		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;}
	};

////
///////
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 );
}

template<class sq_T>
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 ) {
	os << "A:" << ml.A << endl;
	os << "mu:" << ml.mu_const << endl;
	os << "R:" << ml._R() << endl;
	return os;
};

}
#endif //EF_H