/*!
  \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 {
			bdm_error ("Not implemented");
			return 0.0;
		}

		//!Evaluate normalized log-probability
		virtual double evallog (const vec &val) const {
			double tmp;
			tmp = evallog_nn (val) - lognc();
			return tmp;
		}
		//!Evaluate normalized log-probability for many samples
		virtual vec evallog_m (const mat &Val) const {
			vec x (Val.cols());
			for (int i = 0;i < Val.cols();i++) {x (i) = evallog_nn (Val.get_col (i)) ;}
			return x -lognc();
		}
		//!Evaluate normalized log-probability for many samples
		virtual vec evallog_m (const Array<vec> &Val) const {
			vec x (Val.length());
			for (int i = 0;i < Val.length();i++) {x (i) = evallog_nn (Val (i)) ;}
			return x -lognc();
		}

		//!Power of the density, used e.g. to flatten the density
		virtual void pow (double p) {
			bdm_error ("Not implemented");
		}
};


//! Estimator for Exponential family
class BMEF : public BM
{
	protected:
		//! forgetting factor
		double frg;
		//! cached value of lognc() in the previous step (used in evaluation of \c ll )
		double last_lognc;
	public:
		//! Default constructor (=empty constructor)
		BMEF (double frg0 = 1.0) : BM (), frg (frg0) {}
		//! Copy constructor
		BMEF (const BMEF &B) : BM (B), frg (B.frg), last_lognc (B.last_lognc) {}
		//!get statistics from another model
		virtual void set_statistics (const BMEF* BM0) {
			bdm_error ("Not implemented");
		}

		//! Weighted update of sufficient statistics (Bayes rule)
		virtual void bayes (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) {
			bdm_error ("Not implemented");
		}

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

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

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

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

		enorm () : eEF (), mu (), R () {};
		enorm (const vec &mu, const sq_T &R) {set_parameters (mu, R);}
		void set_parameters (const vec &mu, const sq_T &R);
		/*! Create Normal density
		\f[ f(rv) = N(\mu, R) \f]
		from structure
		\code
		class = 'enorm<ldmat>', (OR) 'enorm<chmat>', (OR) 'enorm<fsqmat>';
		mu    = [];                  // mean value
		R     = [];                  // variance, square matrix of appropriate dimension
		\endcode
		*/
		void from_setting (const Setting &root);
		void validate() {
			bdm_assert (mu.length() == R.rows(), "mu and R parameters do not match");
			dim = mu.length();
		}
		//!@}

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

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

		vec sample() const;

		double evallog_nn (const vec &val) const;
		double lognc () const;
		vec mean() const {return mu;}
		vec variance() const {return diag (R.to_mat());}
//	mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
		shared_ptr<mpdf> 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, mpdf &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 );


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

		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;

		//! expected values of the linear coefficient and the covariance matrix are written to \c M and \c R , respectively
		void mean_mat (mat &M, mat&R) const;
		//! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
		double evallog_nn (const vec &val) const;
		double lognc () const;
		void pow (double p) {V *= p;nu *= p;};

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

		ldmat& _V() {return V;}
		const ldmat& _V() const {return V;}
		double& _nu()  {return nu;}
		const double& _nu() const {return nu;}
		const int & _dimx() const {return dimx;}
		/*! Create Gauss-inverse-Wishart density 
		\f[ f(rv) = GiW(V,\nu) \f]
		from structure
		\code
		class = 'egiw';
		V     = [];               // square matrix
		dV    = [];               // vector of diagonal of V (when V not given)
		nu    = [];               // scalar \nu ((almost) degrees of freedom)
								  // when missing, it will be computed to obtain proper pdf
		dimx  = [];               // dimension of the wishart part
		rv = RV({'name'})         // description of RV
		rvc = RV({'name'})        // description of RV in condition
		\endcode
		*/
		
		void from_setting (const Setting &set) {
			epdf::from_setting(set);
			UI::get (dimx, set, "dimx", UI::compulsory);
			if (!UI::get (nu, set, "nu", UI::optional)) {nu=-1;}
			mat V;
			if (!UI::get (V, set, "V", UI::optional)){
				vec dV;
				UI::get (dV, set, "dV", UI::compulsory);
				set_parameters (dimx, ldmat(dV), nu);
				
			} else {
				set_parameters (dimx, V, nu);
			}
		}
		void validate(){
			// check sizes, rvs etc.
		}
		//!@}
};
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);};
		eDirich (const vec &beta0) {set_parameters (beta0);};
		void set_parameters (const vec &beta0) {
			beta = beta0;
			dim = beta.length();
		}
		//!@}

		//! using sampling procedure from wikipedia
		vec sample() const {
			vec y(beta.length());
			for (int i=0; i<beta.length(); i++){
				GamRNG.setup(beta(i),1);
				#pragma omp critical
				y(i)=GamRNG();
			}
			return y/sum(y);
		}

		vec mean() const {return beta / sum (beta);};
		vec variance() const {double gamma = sum (beta); return elem_mult (beta, (gamma-beta)) / (gamma*gamma* (gamma + 1));}
		//! In this instance, val is ...
		double evallog_nn (const vec &val) const {
			double tmp; tmp = (beta - 1) * log (val);
			return tmp;
		}

		double lognc () const {
			double tmp;
			double gam = sum (beta);
			double lgb = 0.0;
			for (int i = 0;i < beta.length();i++) {lgb += lgamma (beta (i));}
			tmp = lgb - lgamma (gam);
			return tmp;
		}

		//!access function
		vec& _beta()  {return beta;}
		/*! configuration structure 
		\code 
		class = 'eDirich';    
		beta  = [];           //parametr beta
		\endcode
		*/
		void from_setting(const Setting &set){
			epdf::from_setting(set);
			UI::get(beta,set, "beta", UI::compulsory);
			validate();
		}
		void validate() {
			//check rv
			dim = beta.length();
		}
};
UIREGISTER(eDirich);

/*! Random Walk on Dirichlet
Using simple assignment 
 \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 mpdf_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(): mpdf_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) {
			mpdf::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() { 
			iepdf.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 &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();}
		}
		//! return correctly typed posterior (covariant return)
		const eDirich& posterior() const {return est;};
		//! constructor function		
		void set_parameters (const vec &beta0) {
			est.set_parameters (beta0);
			if (evalll) {last_lognc = est.lognc();}
		}
};

/*!
 \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;
		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);
			dim = low.length();
		}
		//!@}

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

			UI::get (high, set, "high", UI::compulsory);
			UI::get (low, set, "low", UI::compulsory);
			set_parameters(low,high);
			validate();
		}
		void validate() {
			bdm_assert(high.length()==low.length(), "Incompatible high and low vectors");
			dim = high.length();
			bdm_assert (min (distance) > 0.0, "bad support");
		}
};
UIREGISTER(euni);

//! Uniform density with conditional mean value
class mguni : public mpdf_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){
			mpdf::from_setting(set); //reads rv and rvc
			UI::get(delta,set,"delta",UI::compulsory);
			mean = UI::build<fnc>(set,"mean",UI::compulsory);
			
			iepdf.set_parameters(-delta,delta);
			dimc = mean->dimensionc();
			validate();
		}
};
UIREGISTER(mguni);
/*!
 \brief Normal distributed linear function with linear function of mean value;

 Mean value \f$ \mu=A*\mbox{rvc}+\mu_0 \f$.
*/
template < class sq_T, template <typename> class TEpdf = enorm >
class mlnorm : public mpdf_internal< TEpdf<sq_T> >
{
	protected:
		//! Internal epdf that arise by conditioning on \c rvc
		mat A;
		//! Constant additive term 
		vec mu_const;
//			vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
	public:
		//! \name Constructors
		//!@{
		mlnorm() : mpdf_internal< TEpdf<sq_T> >() {};
		mlnorm (const mat &A, const vec &mu0, const sq_T &R) : mpdf_internal< TEpdf<sq_T> >() {
			set_parameters (A, mu0, R);
		}

		//! Set \c A and \c R
		void set_parameters (const  mat &A0, const vec &mu0, const sq_T &R0) {	
			this->iepdf.set_parameters (zeros (A0.rows()), R0);
			A = A0;
			mu_const = mu0;
			this->dimc = A0.cols();
		}
		//!@}
		//! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
		void condition (const vec &cond) {
			this->iepdf._mu() = A * cond + mu_const;
//R is already assigned;
		}

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

		//! 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) {
			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);
			validate();
		};
		void validate() {
			bdm_assert (A.rows() == mu_const.length(), "mlnorm: A vs. mu mismatch");
			bdm_assert (A.rows() == _R().rows(), "mlnorm: A vs. R mismatch");
			
		}
};
UIREGISTER2 (mlnorm,ldmat);
SHAREDPTR2 ( mlnorm, ldmat );
UIREGISTER2 (mlnorm,fsqmat);
SHAREDPTR2 ( mlnorm, fsqmat );
UIREGISTER2 (mlnorm, chmat);
SHAREDPTR2 ( mlnorm, chmat );

//! Mpdf with general function for mean value
template<class sq_T>
class mgnorm : public mpdf_internal< enorm< sq_T > >
{
	private:
//			vec &mu; WHY NOT?
		shared_ptr<fnc> g;

	public:
		//!default constructor
		mgnorm() : mpdf_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) {
			mpdf::from_setting(set);
			shared_ptr<fnc> g = UI::build<fnc> (set, "g", UI::compulsory);

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

			set_parameters (g, R);
			validate();
		}
		void validate() {
			bdm_assert(g->dimension()==this->dimension(),"incompatible function");
		}
};

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


/*! (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) {
			iepdf._mu() = A * cond + mu_const;
			double zeta;
			//ugly hack!
			if ( (cond.length() + 1) == Lambda.rows()) {
				zeta = Lambda.invqform (concat (cond, vec_1 (1.0)));
			} else {
				zeta = Lambda.invqform (cond);
			}
			_R = Re;
			_R *= (1 + zeta);// / ( nu ); << nu is in Re!!!!!!
		};

		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 mpdf_internal<egamma>
{
	protected:

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

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

	public:
		//! Constructor
		mgamma() : mpdf_internal<egamma>(), k (0),
				_beta (iepdf._beta()) {
		}

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

		void condition (const vec &val) {_beta = k / val;};
		/*! 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) {
			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);
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 mpdf_internal<eigamma>
{
	protected:
		//! Constant \f$k\f$
		double k;

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

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

	public:
		//! \name Constructors
		//!@{
		migamma() : mpdf_internal<eigamma>(),
				k (0),
				_alpha (iepdf._alpha()),
				_beta (iepdf._beta()) {
		}

		migamma (const migamma &m) : mpdf_internal<eigamma>(),
				k (0),
				_alpha (iepdf._alpha()),
				_beta (iepdf._beta()) {
		}
		//!@}

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


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

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

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

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

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


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

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

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

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

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


		/*! 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 mpdf_internal<elognorm>
{
	protected:
		//! parameter 1/2*sigma^2
		double sig2;

		//! access
		vec &mu;
	public:
		//! Constructor
		mlognorm() : mpdf_internal<elognorm>(),
				sig2 (0),
				mu (iepdf._mu()) {
		}

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

			dimc = size;
		};

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

		/*! Create logNormal random Walk
		\f[ f(rv|rvc) = log\mathcal{N}( \log(rvc)-0.5\log(k^2+1), k I) \f]
		from structure
		\code
		class = 'mlognorm';
		k   = 0.1;               // "variance" k
		mu0 = 0.1;               // Initial value of mean
		rv  = RV({'name'})       // description of RV
		rvc = RV({'name'})       // description of RV in condition
		\endcode
		*/
		void from_setting (const Setting &set);

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

};

UIREGISTER (mlognorm);
SHAREDPTR (mlognorm);

/*! inverse Wishart density defined on Choleski decomposition

*/
class eWishartCh : public epdf
{
	protected:
		//! Upper-Triagle of Choleski decomposition of \f$ \Psi \f$
		chmat Y;
		//! dimension of matrix  \f$ \Psi \f$
		int p;
		//! degrees of freedom  \f$ \nu \f$
		double delta;
	public:
		//! Set internal structures
		void set_parameters (const mat &Y0, const double delta0) {Y = chmat (Y0);delta = delta0; p = Y.rows(); dim = p * p; }
		//! 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);
		}
		//! 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);
			dim = W.dimension(); p = Y0.rows();
		}
		vec sample() const {mat iCh; iCh = inv (W.sample_mat()); return vec (iCh._data(), dim);}
		//! access function
		void _setY (const vec &y0) {
			mat Ch (p, p);
			mat iCh (p, p);
			copy_vector (dim, y0._data(), Ch._data());

			iCh = inv (Ch);
			W.setY (iCh);
		}
		virtual double evallog (const vec &val) const {
			chmat X (p);
			const chmat& Y = W.getY();

			copy_vector (p*p, val._data(), X._Ch()._data());
			chmat iX (p);X.inv (iX);
			// compute
// 				\frac{ |\Psi|^{m/2}|X|^{-(m+p+1)/2}e^{-tr(\Psi X^{-1})/2} }{ 2^{mp/2}\Gamma_p(m/2)},
			mat M = Y.to_mat() * iX.to_mat();

			double log1 = 0.5 * p * (2 * Y.logdet()) - 0.5 * (delta + p + 1) * (2 * X.logdet()) - 0.5 * trace (M);
			//Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!

			/*				if (0) {
								mat XX=X.to_mat();
								mat YY=Y.to_mat();

								double log2 = 0.5*p*log(det(YY))-0.5*(delta+p+1)*log(det(XX))-0.5*trace(YY*inv(XX));
								cout << log1 << "," << log2 << endl;
							}*/
			return log1;
		};

};

//! Random Walk on inverse Wishart
class rwiWishartCh : public mpdf_internal<eiWishartCh>
{
	protected:
		//!square root of \f$ \nu-p-1 \f$ - needed for computation of \f$ \Psi \f$ from conditions
		double sqd;
		//!reference point for diagonal
		vec refl;
		//! power of the reference
		double l;
		//! dimension
		int p;

	public:
		rwiWishartCh() : sqd (0), l (0), p (0) {}
		//! constructor function
		void set_parameters (int p0, double k, vec ref0, double l0) {
			p = p0;
			double delta = 2 / (k * k) + p + 3;
			sqd = sqrt (delta - p - 1);
			l = l0;
			refl = pow (ref0, 1 - l);

			iepdf.set_parameters (eye (p), delta);
			dimc = iepdf.dimension();
		}
		void condition (const vec &c) {
			vec z = c;
			int ri = 0;
			for (int i = 0;i < p*p;i += (p + 1)) {//trace diagonal element
				z (i) = pow (z (i), l) * refl (ri);
				ri++;
			}

			iepdf._setY (sqd*z);
		}
};

//! 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);};
		//! Set samples
		void set_parameters (const Array<vec> &Av) {
			bdm_assert(Av.size()>0,"Empty samples"); 
			n = Av.size(); 
			epdf::set_parameters(Av(0).length());
			w=1/n*ones(n);
			samples=Av;
		};
		//! Potentially dangerous, use with care.
		vec& _w()  {return w;};
		//! Potentially dangerous, use with care.
		const vec& _w() const {return w;};
		//! access function
		Array<vec>& _samples() {return samples;};
		//! access function
		const vec& _sample(int i) const {return samples(i);};
		//! access function
		const Array<vec>& _samples() const {return samples;};
		//! Function performs resampling, i.e. removal of low-weight samples and duplication of high-weight samples such that the new samples represent the same density.
		//! The vector with indeces of new samples is returned in variable \c index. 
		void resample ( ivec &index, RESAMPLING_METHOD method = SYSTEMATIC);

		//! Resampling without returning index of new particles. 
		void resample (RESAMPLING_METHOD method = SYSTEMATIC){ivec ind; resample(ind,method);};
		
		//! inherited operation : NOT implemented
		vec sample() const {
			bdm_error ("Not implemented");
			return vec();
		}

		//! inherited operation : NOT implemented
		double evallog (const vec &val) const {
			bdm_error ("Not implemented");
			return 0.0;
		}

		vec mean() const {
			vec pom = zeros (dim);
			for (int i = 0;i < n;i++) {pom += samples (i) * w (i);}
			return pom;
		}
		vec variance() const {
			vec pom = zeros (dim);
			for (int i = 0;i < n;i++) {pom += pow (samples (i), 2) * w (i);}
			return pom -pow (mean(), 2);
		}
		//! For this class, qbounds are minimum and maximum value of the population!
		void qbounds (vec &lb, vec &ub, double perc = 0.95) const {
			// lb in inf so than it will be pushed below;
			lb.set_size (dim);
			ub.set_size (dim);
			lb = std::numeric_limits<double>::infinity();
			ub = -std::numeric_limits<double>::infinity();
			int j;
			for (int i = 0;i < n;i++) {
				for (j = 0;j < dim; j++) {
					if (samples (i) (j) < lb (j)) {lb (j) = samples (i) (j);}
					if (samples (i) (j) > ub (j)) {ub (j) = samples (i) (j);}
				}
			}
		}
};


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

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>
// 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<mpdf> enorm<sq_T>::condition ( const RV &rvn ) const
{
	mlnorm<sq_T> *tmp = new mlnorm<sq_T> ();
	shared_ptr<mpdf> narrow(tmp);
	condition ( rvn, *tmp );
	return narrow;
}

template<class sq_T>
void enorm<sq_T>::condition ( const RV &rvn, mpdf &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);
}

////
///////
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