/*!
  \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;
			//!TODO decide if it is really needed
			virtual void dupdate ( mat &v ) {it_error ( "Not implemented" );};
			//!Evaluate normalized log-probability
			virtual double evallog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
			//!Evaluate normalized log-probability
			virtual double evallog ( const vec &val ) const {
				double tmp;
				tmp= evallog_nn ( val )-lognc();
// 				it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); 
				return tmp;}
			//!Evaluate normalized log-probability for many samples
			virtual vec evallog ( 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();
			}
			//!Power of the density, used e.g. to flatten the density
			virtual void pow ( double p ) {it_error ( "Not implemented" );};
	};

	/*!
	* \brief Exponential family model.

	* More?...
	*/

	class mEF : public mpdf
	{

		public:
			//! Default constructor
			mEF ( ) :mpdf ( ) {};
	};

//! Estimator for Exponential family
	class BMEF : public BM
	{
		protected:
			//! forgetting factor
			double frg;
			//! cached value of lognc() in the previous step (used in evaluation of \c ll )
			double last_lognc;
		public:
			//! Default constructor (=empty constructor)
			BMEF ( double frg0=1.0 ) :BM ( ), frg ( frg0 ) {}
			//! Copy constructor
			BMEF ( const BMEF &B ) :BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
			//!get statistics from another model
			virtual void set_statistics ( const BMEF* BM0 ) {it_error ( "Not implemented" );};
			//! Weighted update of sufficient statistics (Bayes rule)
			virtual void bayes ( const vec &data, const double w ) {};
			//original Bayes
			void bayes ( const vec &dt );
			//!Flatten the posterior according to the given BMEF (of the same type!)
			virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
			//!Flatten the posterior as if to keep nu0 data
//	virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}

			BMEF* _copy_ () const {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
	};

	template<class sq_T>
	class mlnorm;

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

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

			enorm ( ) :eEF ( ), mu ( ),R ( ) {};
			enorm ( const vec &mu,const sq_T &R ) {set_parameters ( mu,R );}
			void set_parameters ( const vec &mu,const sq_T &R );
			void from_setting(const Setting &root);
			void validate() {
				it_assert(mu.length()==R.rows(),"parameters mismatch");
				dim = mu.length();
			}
			//!@}

			//! \name Mathematical operations
			//!@{

			//! dupdate in exponential form (not really handy)
			void dupdate ( mat &v,double nu=1.0 );

			vec sample() const;

			double evallog_nn ( const vec &val ) const;
			double lognc () const;
			vec mean() const {return mu;}
			vec variance() const {return diag ( R.to_mat() );}
//	mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
			mpdf* condition ( const RV &rvn ) const ;
			enorm<sq_T>* marginal ( const RV &rv ) const;
//			epdf* marginal ( const RV &rv ) const;
			//!@}

			//! \name Access to attributes
			//!@{

			vec& _mu() {return mu;}
			void set_mu ( const vec mu0 ) { mu=mu0;}
			sq_T& _R() {return R;}
			const sq_T& _R() const {return R;}
			//!@}

	};
	UIREGISTER(enorm<chmat>);
	UIREGISTER(enorm<ldmat>);
	UIREGISTER(enorm<fsqmat>);


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

	* For \f$p\f$-variate densities, given rv.count() should be \f$p\times\f$ V.rows().
	*
	*/
	class egiw : public eEF
	{
		protected:
			//! Extended information matrix of sufficient statistics
			ldmat V;
			//! Number of data records (degrees of freedom) of sufficient statistics
			double nu;
			//! Dimension of the output
			int dimx;
			//! Dimension of the regressor
			int nPsi;
		public:
			//!\name Constructors
			//!@{
			egiw() :eEF() {};
			egiw ( int dimx0, ldmat V0, double nu0=-1.0 ) :eEF() {set_parameters ( dimx0,V0, nu0 );};

			void set_parameters ( int dimx0, ldmat V0, double nu0=-1.0 )
			{
				dimx=dimx0;
				nPsi = V0.rows()-dimx;
				dim = dimx* ( dimx+nPsi ); // size(R) + size(Theta)

				V=V0;
				if ( nu0<0 )
				{
					nu = 0.1 +nPsi +2*dimx +2; // +2 assures finite expected value of R
					// terms before that are sufficient for finite normalization
				}
				else
				{
					nu=nu0;
				}
			}
			//!@}

			vec sample() const;
			vec mean() const;
			vec variance() const;

			//! LS estimate of \f$\theta\f$
			vec est_theta() const;

			//! Covariance of the LS estimate
			ldmat est_theta_cov() const;

			void mean_mat ( mat &M, mat&R ) const;
			//! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
			double evallog_nn ( const vec &val ) const;
			double lognc () const;
			void pow ( double p ) {V*=p;nu*=p;};

			//! \name Access attributes
			//!@{

			ldmat& _V() {return V;}
			const ldmat& _V() const {return V;}
			double& _nu()  {return nu;}
			const double& _nu() const {return nu;}
			void from_setting(const Setting &set){
				UI::get(nu, set, "nu", UI::compulsory);
				UI::get(dimx, set, "dimx", UI::compulsory);
				mat V;
				UI::get(V,set,"V", UI::compulsory);
				set_parameters(dimx, V, nu);
				RV* rv=UI::build<RV>(set,"rv", UI::compulsory);
				set_rv(*rv);
				delete rv;
			}
			//!@}
	};
	UIREGISTER(egiw);

	/*! \brief Dirichlet posterior density

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

			eDirich () : eEF ( ) {};
			eDirich ( const eDirich &D0 ) : eEF () {set_parameters ( D0.beta );};
			eDirich ( const vec &beta0 ) {set_parameters ( beta0 );};
			void set_parameters ( const vec &beta0 )
			{
				beta= beta0;
				dim = beta.length();
			}
			//!@}

			vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
			vec mean() const {return beta/sum(beta);};
			vec variance() const {double gamma =sum(beta); return elem_mult ( beta, ( beta+1 ) ) / ( gamma* ( gamma+1 ) );}
			//! In this instance, val is ...
			double evallog_nn ( const vec &val ) const
			{
				double tmp; tmp= ( beta-1 ) *log ( val );
// 				it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
				return tmp;
			};
			double lognc () const
			{
				double tmp;
				double gam=sum ( beta );
				double lgb=0.0;
				for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
				tmp= lgb-lgamma ( gam );
// 				it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
				return tmp;
			};
			//!access function
			vec& _beta()  {return beta;}
			//!Set internal parameters
	};

//! \brief Estimator for Multinomial density
	class multiBM : public BMEF
	{
		protected:
			//! Conjugate prior and posterior
			eDirich est;
			//! Pointer inside est to sufficient statistics
			vec &beta;
		public:
			//!Default constructor
			multiBM ( ) : BMEF ( ),est ( ),beta ( est._beta() )
			{
				if ( beta.length() >0 ) {last_lognc=est.lognc();}
				else{last_lognc=0.0;}
			}
			//!Copy constructor
			multiBM ( const multiBM &B ) : BMEF ( B ),est ( B.est ),beta ( est._beta() ) {}
			//! Sets sufficient statistics to match that of givefrom mB0
			void set_statistics ( const BM* mB0 ) {const multiBM* mB=dynamic_cast<const multiBM*> ( mB0 ); beta=mB->beta;}
			void bayes ( const vec &dt )
			{
				if ( frg<1.0 ) {beta*=frg;last_lognc=est.lognc();}
				beta+=dt;
				if ( evalll ) {ll=est.lognc()-last_lognc;}
			}
			double logpred ( const vec &dt ) const
			{
				eDirich pred ( est );
				vec &beta = pred._beta();

				double lll;
				if ( frg<1.0 )
					{beta*=frg;lll=pred.lognc();}
				else
					if ( evalll ) {lll=last_lognc;}
					else{lll=pred.lognc();}

				beta+=dt;
				return pred.lognc()-lll;
			}
			void flatten ( const BMEF* B )
			{
				const multiBM* E=dynamic_cast<const multiBM*> ( B );
				// sum(beta) should be equal to sum(B.beta)
				const vec &Eb=E->beta;//const_cast<multiBM*> ( E )->_beta();
				beta*= ( sum ( Eb ) /sum ( beta ) );
				if ( evalll ) {last_lognc=est.lognc();}
			}
			const epdf& posterior() const {return est;};
			const eDirich* _e() const {return &est;};
			void set_parameters ( const vec &beta0 )
			{
				est.set_parameters ( beta0 );
				if ( evalll ) {last_lognc=est.lognc();}
			}
	};

	/*!
	 \brief Gamma posterior density

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

	class egamma : public eEF
	{
		protected:
			//! Vector \f$\alpha\f$
			vec alpha;
			//! Vector \f$\beta\f$
			vec beta;
		public :
			//! \name Constructors
			//!@{
			egamma ( ) :eEF ( ), alpha ( 0 ), beta ( 0 ) {};
			egamma ( const vec &a, const vec &b ) {set_parameters ( a, b );};
			void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;dim = alpha.length();};
			//!@}

			vec sample() const;
			//! TODO: is it used anywhere?
//	mat sample ( int N ) const;
			double evallog ( const vec &val ) const;
			double lognc () const;
			//! Returns poiter to alpha and beta. Potentially dengerous: use with care!
			vec& _alpha() {return alpha;}
			vec& _beta() {return beta;}
			vec mean() const {return elem_div ( alpha,beta );}
			vec variance() const {return elem_div ( alpha,elem_mult ( beta,beta ) ); }
			
			//! Load from structure with elements:
			//!  \code
			//! { alpha = [...];         // vector of alpha
			//!   beta = [...];          // vector of beta
			//!   rv = {class="RV",...}  // description
			//! }
			//! \endcode
			//!@}
			void from_setting(const Setting &set){
				epdf::from_setting(set); // reads rv
				UI::get(alpha,set,"alpha", UI::compulsory);
				UI::get(beta,set,"beta", UI::compulsory);
				validate();
			}
			void validate(){
				it_assert(alpha.length() ==beta.length(), "parameters do not match");
				dim =alpha.length();
			}
	};
UIREGISTER(egamma);
	/*!
	 \brief Inverse-Gamma posterior density

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

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

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

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

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

//!  Uniform distributed density on a rectangular support

	class euni: public epdf
	{
		protected:
//! lower bound on support
			vec low;
//! upper bound on support
			vec high;
//! internal
			vec distance;
//! normalizing coefficients
			double nk;
//! cache of log( \c nk )
			double lnk;
		public:
			//! \name Constructors
			//!@{
			euni ( ) :epdf ( ) {}
			euni ( const vec &low0, const vec &high0 ) {set_parameters ( low0,high0 );}
			void set_parameters ( const vec &low0, const vec &high0 )
			{
				distance = high0-low0;
				it_assert_debug ( min ( distance ) >0.0,"bad support" );
				low = low0;
				high = high0;
				nk = prod ( 1.0/distance );
				lnk = log ( nk );
				dim = low.length();
			}
			//!@}

			double eval ( const vec &val ) const  {return nk;}
			double evallog ( const vec &val ) const  {
				if (any(val<low) && any(val>high)) {return inf;}
				else return lnk;
			}
			vec sample() const
			{
				vec smp ( dim );
#pragma omp critical
				UniRNG.sample_vector ( dim ,smp );
				return low+elem_mult ( distance,smp );
			}
			//! set values of \c low and \c high
			vec mean() const {return ( high-low ) /2.0;}
			vec variance() const {return ( pow ( high,2 ) +pow ( low,2 ) +elem_mult ( high,low ) ) /3.0;}
			//! Load from structure with elements:
			//!  \code
			//! { high = [...];          // vector of upper bounds
			//!   low = [...];           // vector of lower bounds
			//!   rv = {class="RV",...}  // description of RV
			//! }
			//! \endcode
			//!@}
			void from_setting(const Setting &set){
				epdf::from_setting(set); // reads rv and rvc

				UI::get(high,set,"high", UI::compulsory);
				UI::get(low,set,"low", UI::compulsory);
			}
	};


	/*!
	 \brief Normal distributed linear function with linear function of mean value;

	 Mean value \f$mu=A*rvc+mu_0\f$.
	*/
	template<class sq_T>
	class mlnorm : public mEF
	{
		protected:
			//! Internal epdf that arise by conditioning on \c rvc
			shared_ptr<enorm<sq_T> > iepdf;
			mat A;
			vec mu_const;
			vec& _mu; //cached epdf.mu;
		public:
			//! \name Constructors
			//!@{
			mlnorm():iepdf(new enorm<sq_T>()), _mu(iepdf->_mu()) { set_ep(iepdf); };
			mlnorm (const mat &A, const vec &mu0, const sq_T &R ) :iepdf(new enorm<sq_T>()), _mu(iepdf->_mu())
			{
				set_ep(iepdf); set_parameters ( A,mu0,R );
			}

			//! Set \c A and \c R
			void set_parameters ( const  mat &A, const vec &mu0, const sq_T &R );
			//!@}
			//! 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 );

			//!access function
			vec& _mu_const() {return mu_const;}
			//!access function
			mat& _A() {return A;}
			//!access function
			mat _R() { return iepdf->_R().to_mat(); }

			template<class sq_M>
			friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M> &ml );
			
			void from_setting(const Setting &set){
				mpdf::from_setting(set);	

				UI::get(A,set,"A", UI::compulsory);
				UI::get(mu_const,set,"const", UI::compulsory);
				mat R0;
				UI::get(R0,set,"R", UI::compulsory);
				set_parameters(A,mu_const,R0);
			};
	};
	UIREGISTER(mlnorm<ldmat>);
	UIREGISTER(mlnorm<fsqmat>);
	UIREGISTER(mlnorm<chmat>);

//! Mpdf with general function for mean value
	template<class sq_T>
	class mgnorm : public mEF
	{
		protected:
			//! Internal epdf that arise by conditioning on \c rvc
			shared_ptr<enorm<sq_T> > iepdf;
			vec &mu;
			fnc* g;
		public:
			//!default constructor
			mgnorm():iepdf(new enorm<sq_T>()), mu(iepdf->_mu()) { set_ep(iepdf); }
			//!set mean function
			void set_parameters ( fnc* g0, const sq_T &R0 ) {g=g0; iepdf->set_parameters ( zeros ( g->dimension() ), R0 );}
			void condition ( const vec &cond ) {mu=g->eval ( cond );};


			/*! UI for mgnorm

			The mgnorm is constructed from a structure with fields:
			\code
			system = {
				type = "mgnorm";
				// function for mean value evolution
				g = {type="fnc"; ... }

				// variance
				R = [1, 0,
					 0, 1];
				// --OR --
				dR = [1, 1];

				// == OPTIONAL ==

				// description of y variables
				y = {type="rv"; names=["y", "u"];};
				// description of u variable
				u = {type="rv"; names=[];}
			};
			\endcode

			Result if
			*/

			void from_setting( const Setting &set ) 
			{	
				fnc* g = UI::build<fnc>( set, "g", UI::compulsory );

				mat R;					
				vec dR;					
				if ( UI::get( dR, set, "dR" ) )
					R=diag(dR);
				else 
					UI::get( R, set, "R", UI::compulsory); 
		
				set_parameters(g,R);
			}
	};

	UIREGISTER(mgnorm<chmat>);


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

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

	Perhaps a moment-matching technique?
	*/
	class mlstudent : public mlnorm<ldmat>
	{
		protected:
			ldmat Lambda;
			ldmat &_R;
			ldmat Re;
		public:
			mlstudent ( ) :mlnorm<ldmat> (),
					Lambda (),	_R ( iepdf->_R() ) {}
			void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 )
			{
				it_assert_debug ( A0.rows() ==mu0.length(),"" );
				it_assert_debug ( R0.rows() ==A0.rows(),"" );

				iepdf->set_parameters ( mu0,Lambda ); //
				A = A0;
				mu_const = mu0;
				Re=R0;
				Lambda = Lambda0;
			}
			void condition ( const vec &cond )
			{
				_mu = A*cond + mu_const;
				double zeta;
				//ugly hack!
				if ( ( cond.length() +1 ) ==Lambda.rows() )
				{
					zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) );
				}
				else
				{
					zeta = Lambda.invqform ( cond );
				}
				_R = Re;
				_R*= ( 1+zeta );// / ( nu ); << nu is in Re!!!!!!
			};

	};
	/*!
	\brief  Gamma random walk

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

	The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
	*/
	class mgamma : public mEF
	{
		protected:
			//! Internal epdf that arise by conditioning on \c rvc
			shared_ptr<egamma> iepdf;

			//! Constant \f$k\f$
			double k;

			//! cache of iepdf.beta
			vec &_beta;

		public:
			//! Constructor
			mgamma():iepdf(new egamma()), k(0),
			    _beta(iepdf->_beta()) {
			    set_ep(iepdf);
			}

			//! Set value of \c k
			void set_parameters(double k, const vec &beta0);

			void condition ( const vec &val ) {_beta=k/val;};
			//! Load from structure with elements:
			//!  \code
			//! { alpha = [...];         // vector of alpha
			//!   k = 1.1;               // multiplicative constant k
			//!   rv = {class="RV",...}  // description of RV
			//!   rvc = {class="RV",...} // description of RV in condition
			//! }
			//! \endcode
			//!@}
			void from_setting(const Setting &set){
				mpdf::from_setting(set); // reads rv and rvc
				vec betatmp; // ugly but necessary
				UI::get(betatmp,set,"beta", UI::compulsory);
				UI::get(k,set,"k", UI::compulsory);
				set_parameters(k,betatmp);
			}
	};
	UIREGISTER(mgamma);
	
	/*!
	\brief  Inverse-Gamma random walk

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

	The standard deviation of the walk is then: \f$ \mu/\sqrt(k)\f$.
	 */
	class migamma : public mEF
	{
		protected:
			//! Internal epdf that arise by conditioning on \c rvc
			shared_ptr<eigamma> iepdf;

			//! Constant \f$k\f$
			double k;

			//! cache of iepdf.alpha
			vec &_alpha;

			//! cache of iepdf.beta
			vec &_beta;

		public:
			//! \name Constructors
			//!@{
			migamma():iepdf(new eigamma()),
			    k(0),
			    _alpha(iepdf->_alpha()),
			    _beta(iepdf->_beta()) {
			    set_ep(iepdf);
			}

			migamma(const migamma &m):iepdf(m.iepdf),
			    k(0),
			    _alpha(iepdf->_alpha()),
			    _beta(iepdf->_beta()) {
			    set_ep(iepdf);
			}
			//!@}

			//! Set value of \c k
			void set_parameters ( int len, double k0 )
			{
				k=k0;
				iepdf->set_parameters ( ( 1.0/ ( k*k ) +2.0 ) *ones ( len ) /*alpha*/, ones ( len ) /*beta*/ );
				dimc = dimension();
			};
			void condition ( const vec &val )
			{
				_beta=elem_mult ( val, ( _alpha-1.0 ) );
			};
	};


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

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

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

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

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


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

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

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

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

			void condition ( const vec &val )
			{
				vec mean=elem_mult ( refl,pow ( val,l ) );
				migamma::condition ( mean );
			};

			/*! UI for migamma_ref

			The migamma_ref is constructed from a structure with fields:
			\code
			system = {
				type = "migamma_ref";
				ref = [1e-5; 1e-5; 1e-2 1e-3];            // reference vector
				l = 0.999;                                // constant l
				k = 0.1;                                  // constant k
				
				// == OPTIONAL ==
				// description of y variables
				y = {type="rv"; names=["y", "u"];};
				// description of u variable
				u = {type="rv"; names=[];}
			};
			\endcode

			Result if
			 */
			void from_setting( const Setting &set );

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


	UIREGISTER(migamma_ref);

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

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

	*/
	class elognorm: public enorm<ldmat>
	{
		public:
			vec sample() const {return exp ( enorm<ldmat>::sample() );};
			vec mean() const {vec var=enorm<ldmat>::variance();return exp ( mu - 0.5*var );};

	};

	/*!
	\brief  Log-Normal random walk

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

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

	The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
	 */
	class mlognorm : public mpdf
	{
		protected:
			shared_ptr<elognorm> eno;

			//! parameter 1/2*sigma^2
			double sig2;

			//! access
			vec &mu;
		public:
			//! Constructor
			mlognorm():eno(new elognorm()),
			    sig2(0),
			    mu(eno->_mu()) {
			    set_ep(eno);
			}

			//! Set value of \c k
			void set_parameters ( int size, double k )
			{
				sig2 = 0.5*log ( k*k+1 );
				eno->set_parameters ( zeros ( size ),2*sig2*eye ( size ) );

				dimc = size;
			};

			void condition ( const vec &val )
			{
				mu=log ( val )-sig2;//elem_mult ( refl,pow ( val,l ) );
			};

			/*! UI for mlognorm

			The mlognorm is constructed from a structure with fields:
			\code
			system = {
				type = "mlognorm";
				k = 0.1;                                  // constant k
				mu0 = [1., 1.];
				
				// == OPTIONAL ==
				// description of y variables
				y = {type="rv"; names=["y", "u"];};
				// description of u variable
				u = {type="rv"; names=[];}
			};
			\endcode

			 */
			void from_setting( const Setting &set );

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

	};
	
	UIREGISTER(mlognorm);

	/*! inverse Wishart density defined on Choleski decomposition

	*/
	class eWishartCh : public epdf
	{
		protected:
			//! Upper-Triagle of Choleski decomposition of \f$ \Psi \f$
			chmat Y;
			//! dimension of matrix  \f$ \Psi \f$
			int p;
			//! degrees of freedom  \f$ \nu \f$
			double delta;
		public:
			void set_parameters ( const mat &Y0, const double delta0 ) {Y=chmat ( Y0 );delta=delta0; p=Y.rows(); dim = p*p; }
			mat sample_mat() const
			{
				mat X=zeros ( p,p );

				//sample diagonal
				for ( int i=0;i<p;i++ )
				{
					GamRNG.setup ( 0.5* ( delta-i ) , 0.5 ); // no +1 !! index if from 0
#pragma omp critical
					X ( i,i ) =sqrt ( GamRNG() );
				}
				//do the rest
				for ( int i=0;i<p;i++ )
				{
					for ( int j=i+1;j<p;j++ )
					{
#pragma omp critical
						X ( i,j ) =NorRNG.sample();
					}
				}
				return X*Y._Ch();// return upper triangular part of the decomposition
			}
			vec sample () const
			{
				return vec ( sample_mat()._data(),p*p );
			}
			//! fast access function y0 will be copied into Y.Ch.
			void setY ( const mat &Ch0 ) {copy_vector ( dim,Ch0._data(), Y._Ch()._data() );}
			//! fast access function y0 will be copied into Y.Ch.
			void _setY ( const vec &ch0 ) {copy_vector ( dim, ch0._data(), Y._Ch()._data() ); }
			//! access function
			const chmat& getY()const {return Y;}
	};

	class eiWishartCh: public epdf
	{
		protected:
			eWishartCh W;
			int p;
			double delta;
		public:
			void set_parameters ( const mat &Y0, const double delta0) {
				delta = delta0;
				W.set_parameters ( inv ( Y0 ),delta0 ); 
				dim = W.dimension(); p=Y0.rows();
			}
			vec sample() const {mat iCh; iCh=inv ( W.sample_mat() ); return vec ( iCh._data(),dim );}
			void _setY ( const vec &y0 )
			{
				mat Ch ( p,p );
				mat iCh ( p,p );
				copy_vector ( dim, y0._data(), Ch._data() );
				
				iCh=inv ( Ch );
				W.setY ( iCh );
			}
			virtual double evallog ( const vec &val ) const {
				chmat X(p);
				const chmat& Y=W.getY();
				 
				copy_vector(p*p,val._data(),X._Ch()._data());
				chmat iX(p);X.inv(iX);
				// compute  
// 				\frac{ |\Psi|^{m/2}|X|^{-(m+p+1)/2}e^{-tr(\Psi X^{-1})/2} }{ 2^{mp/2}\Gamma_p(m/2)},
				mat M=Y.to_mat()*iX.to_mat();
				
				double log1 = 0.5*p*(2*Y.logdet())-0.5*(delta+p+1)*(2*X.logdet())-0.5*trace(M); 
				//Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!
				
/*				if (0) {
					mat XX=X.to_mat();
					mat YY=Y.to_mat();
					
					double log2 = 0.5*p*log(det(YY))-0.5*(delta+p+1)*log(det(XX))-0.5*trace(YY*inv(XX)); 
					cout << log1 << "," << log2 << endl;
				}*/
				return log1;				
			};
			
	};

	class rwiWishartCh : public mpdf
	{
		protected:
			shared_ptr<eiWishartCh> eiW;
			//!square root of \f$ \nu-p-1 \f$ - needed for computation of \f$ \Psi \f$ from conditions
			double sqd;
			//reference point for diagonal
			vec refl;
			double l;
			int p;

		public:
			rwiWishartCh():eiW(new eiWishartCh()),
			    sqd(0), l(0), p(0) {
			    set_ep(eiW);
			}

			void set_parameters ( int p0, double k, vec ref0, double l0  )
			{
				p=p0;
				double delta = 2/(k*k)+p+3;
				sqd=sqrt ( delta-p-1 );
				l=l0;
				refl=pow(ref0,1-l);
				
				eiW->set_parameters ( eye ( p ),delta );
				dimc=eiW->dimension();
			}
			void condition ( const vec &c ) {
				vec z=c;
				int ri=0;
				for(int i=0;i<p*p;i+=(p+1)){//trace diagonal element
					z(i) = pow(z(i),l)*refl(ri);
					ri++;
				}

				eiW->_setY ( sqd*z );
			}
	};

//! Switch between various resampling methods.
	enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
	/*!
	\brief Weighted empirical density

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

			//! Set samples and weights
			void set_statistics ( const vec &w0, const epdf* pdf0 );
			//! Set samples and weights
			void set_statistics ( const epdf* pdf0 , int n ) {set_statistics ( ones ( n ) /n,pdf0 );};
			//! Set sample
			void set_samples ( const epdf* pdf0 );
			//! Set sample
			void set_parameters ( int n0, bool copy=true ) {n=n0; w.set_size ( n0,copy );samples.set_size ( n0,copy );};
			//! Potentially dangerous, use with care.
			vec& _w()  {return w;};
			//! Potentially dangerous, use with care.
			const vec& _w() const {return w;};
			//! access function
			Array<vec>& _samples() {return samples;};
			//! access function
			const Array<vec>& _samples() const {return samples;};
			//! Function performs resampling, i.e. removal of low-weight samples and duplication of high-weight samples such that the new samples represent the same density.
			ivec resample ( RESAMPLING_METHOD method=SYSTEMATIC );
			//! inherited operation : NOT implemneted
			vec sample() const {it_error ( "Not implemented" );return 0;}
			//! inherited operation : NOT implemneted
			double evallog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
			vec mean() const
			{
				vec pom=zeros ( dim );
				for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
				return pom;
			}
			vec variance() const
			{
				vec pom=zeros ( dim );
				for ( int i=0;i<n;i++ ) {pom+=pow ( samples ( i ),2 ) *w ( i );}
				return pom-pow ( mean(),2 );
			}
			//! For this class, qbounds are minimum and maximum value of the population!
			void qbounds ( vec &lb, vec &ub, double perc=0.95 ) const
			{
				// lb in inf so than it will be pushed below;
				lb.set_size ( dim );
				ub.set_size ( dim );
				lb = std::numeric_limits<double>::infinity();
				ub = -std::numeric_limits<double>::infinity();
				int j;
				for ( int i=0;i<n;i++ )
				{
					for ( j=0;j<dim; j++ )
					{
						if ( samples ( i ) ( j ) <lb ( j ) ) {lb ( j ) =samples ( i ) ( j );}
						if ( samples ( i ) ( j ) >ub ( j ) ) {ub ( j ) =samples ( i ) ( j );}
					}
				}
			}
	};


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

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

	template<class sq_T>
	void enorm<sq_T>::from_setting(const Setting &set){
		epdf::from_setting(set); //reads rv
		
		UI::get(mu,set,"mu", UI::compulsory);
		mat Rtmp;// necessary for conversion
		UI::get(Rtmp,set,"R", UI::compulsory);
		R=Rtmp; // conversion
		validate();
	}

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

// 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>
	void mlnorm<sq_T>::set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 )
	{
		it_assert_debug ( A0.rows() ==mu0.length(),"" );
		it_assert_debug ( A0.rows() ==R0.rows(),"" );

		iepdf->set_parameters(zeros(A0.rows()), R0);
		A = A0;
		mu_const = mu0;
		dimc = A0.cols();
	}

// 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>
	void mlnorm<sq_T>::condition ( const vec &cond )
	{
		_mu = A*cond + mu_const;
//R is already assigned;
	}

	template<class sq_T>
	enorm<sq_T>* enorm<sq_T>::marginal ( const RV &rvn ) const
	{
		it_assert_debug ( isnamed(), "rv description is not assigned" );
		ivec irvn = rvn.dataind ( rv );

		sq_T Rn ( R,irvn ); //select rows and columns of R

		enorm<sq_T>* tmp = new enorm<sq_T>;
		tmp->set_rv ( rvn );
		tmp->set_parameters ( mu ( irvn ), Rn );
		return tmp;
	}

	template<class sq_T>
	mpdf* enorm<sq_T>::condition ( const RV &rvn ) const
	{

		it_assert_debug ( isnamed(),"rvs are not assigned" );

		RV rvc = rv.subt ( rvn );
		it_assert_debug ( ( rvc._dsize() +rvn._dsize() ==rv._dsize() ),"wrong rvn" );
		//Permutation vector of the new R
		ivec irvn = rvn.dataind ( rv );
		ivec irvc = rvc.dataind ( rv );
		ivec perm=concat ( irvn , irvc );
		sq_T Rn ( R,perm );

		//fixme - could this be done in general for all sq_T?
		mat S=Rn.to_mat();
		//fixme
		int n=rvn._dsize()-1;
		int end=R.rows()-1;
		mat S11 = S.get ( 0,n, 0, n );
		mat S12 = S.get ( 0, n , rvn._dsize(), end );
		mat S22 = S.get ( rvn._dsize(), end, rvn._dsize(), end );

		vec mu1 = mu ( irvn );
		vec mu2 = mu ( irvc );
		mat A=S12*inv ( S22 );
		sq_T R_n ( S11 - A *S12.T() );

		mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( );
		tmp->set_rv ( rvn ); tmp->set_rvc ( rvc );
		tmp->set_parameters ( A,mu1-A*mu2,R_n );
		return tmp;
	}

///////////

	template<class sq_T>
	std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_T> &ml )
	{
		os << "A:"<< ml.A<<endl;
		os << "mu:"<< ml.mu_const<<endl;
		os << "R:" << ml.iepdf->_R().to_mat() <<endl;
		return os;
	};

}
#endif //EF_H