/*!
  \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 "libBM.h"
#include "../math/chmat.h"
//#include <std>

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_ ( bool changerv=false ) {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 );
	//!@}

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

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

	vec sample() const;
	mat sample ( int N ) 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 ;
	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;}
	//!@}

};

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

/*! \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;
	//!speedup variable
	double gamma;
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();
		gamma = sum ( beta );
	}
	//!@}

	vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
	vec mean() const {return beta/gamma;};
	vec variance() const {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& _epdf() 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(), beta() {};
	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 ) ); }
};

/*!
 \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]

 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 eEF {
protected:
	//!internal egamma
	egamma eg;
	//! Vector \f$\alpha\f$
	vec &alpha;
	//! Vector \f$\beta\f$ (in fact it is 1/beta as used in definition of iG)
	vec &beta;
public :
	//! \name Constructors
	//!@{
	eigamma ( ) :eEF ( ), eg(),alpha ( eg._alpha() ), beta ( eg._beta() ) {};
	eigamma ( const vec &a, const vec &b ):eEF ( ), eg(),alpha ( eg._alpha() ), beta ( eg._beta() ) {eg.set_parameters ( a,b );};
	void set_parameters ( const vec &a, const vec &b ) {eg.set_parameters ( a,b );};
	//!@}

	vec sample() const {return 1.0/eg.sample();};
	//! TODO: is it used anywhere?
//	mat sample ( int N ) const;
	double evallog ( const vec &val ) const {return eg.evallog ( val );};
	double lognc () const {return eg.lognc();};
	//! 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 );}
	vec& _alpha() {return alpha;}
	vec& _beta() {return beta;}
};
/*
//! 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  {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;}
};


/*!
 \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
	enorm<sq_T> epdf;
	mat A;
	vec mu_const;
	vec& _mu; //cached epdf.mu;
public:
	//! \name Constructors
	//!@{
	mlnorm ( ) :mEF (),epdf ( ),A ( ),_mu ( epdf._mu() ) {ep =&epdf; };
	mlnorm ( const  mat &A, const vec &mu0, const sq_T &R ) :epdf ( ),_mu ( epdf._mu() ) {
		ep =&epdf; 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 epdf._R().to_mat();}

	template<class sq_M>
	friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M> &ml );
};

/*! (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 ( epdf._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(),"" );

		epdf.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
	egamma epdf;
	//! Constant \f$k\f$
	double k;
	//! cache of epdf.beta
	vec &_beta;

public:
	//! Constructor
	mgamma ( ) : mEF ( ), epdf (), _beta ( epdf._beta() ) {ep=&epdf;};
	//! Set value of \c k
	void set_parameters ( double k, const vec &beta0 );
	void condition ( const vec &val ) {_beta=k/val;};
};

/*!
\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
	eigamma epdf;
	//! Constant \f$k\f$
	double k;
	//! cache of epdf.alpha
	vec &_alpha;
	//! cache of epdf.beta
	vec &_beta;

public:
	//! \name Constructors
	//!@{
	migamma ( ) : mEF (), epdf ( ), _alpha ( epdf._alpha() ), _beta ( epdf._beta() ) {ep=&epdf;};
	migamma ( const migamma &m ) : mEF (), epdf ( m.epdf ), _alpha ( epdf._alpha() ), _beta ( epdf._beta() ) {ep=&epdf;};
	//!@}

	//! Set value of \c k
	void set_parameters ( int len, double k0 ) {
		k=k0;
		epdf.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_fix : public migamma {
protected:
	//! parameter l
	double l;
	//! reference vector
	vec refl;
public:
	//! Constructor
	migamma_fix ( ) : 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 );
	};
};
//! 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 ( ) {};
	eEmp (const eEmp &e ): epdf(e), w(e.w), samples(e.samples) {};
	//!@}
	
	//! Set samples and weights
	void set_parameters ( const vec &w0, const epdf* pdf0 );
	//! Set sample
	void set_samples ( const epdf* pdf0 );
	//! Set sample
	void set_n ( int n0, bool copy=true ) {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 );
	}
};


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

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;
	dim = mu0.length();
};

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>
mat enorm<sq_T>::sample ( int N ) const {
	mat X ( dim,N );
	vec x ( dim );
	vec pom;
	int i;

	for ( i=0;i<N;i++ ) {
#pragma omp critical
		NorRNG.sample_vector ( dim,x );
		pom = R.sqrt_mult ( x );
		pom +=mu;
		X.set_col ( i, pom );
	}

	return X;
};

// 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(),"" );

	epdf.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>
epdf* 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.epdf._R().to_mat() <<endl;
	return os;
};

}
#endif //EF_H
