/*!
  \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 <itpp/itbase.h>
#include "libDC.h"
#include "libBM.h"
//#include <std>

using namespace itpp;

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

* More?...
*/
class eEF : public epdf {

public:
	virtual void tupdate( double phi, mat &vbar, double nubar ) {};
	virtual void dupdate( mat &v,double nu=1.0 ) {};
};

class mEF : public mpdf {

public:

};

/*!
* \brief General exponential family density

* More?...
*/
template<class sq_T>
class enorm : public eEF {
	int dim;
	vec mu;
	sq_T R;
public:
	enorm( RV &rv, vec &mu, sq_T &R );
	enorm();
	void tupdate( double phi, mat &vbar, double nubar );
	void dupdate( mat &v,double nu=1.0 );
	vec sample();
	mat sample(int N);
	double eval( const vec &val );
	Normal_RNG RNG;
};

/*!
 \brief
*/
template<class sq_T>
class mlnorm : public mEF {
	enorm<sq_T> epdf;
	mat A;
public:
	//! Constructor
	mlnorm( RV &rv,RV &rvc, mat &A, sq_T &R );
	//!Generate one sample of the posterior
	vec samplecond( vec &cond, double &lik );
	mat samplecond( vec &cond, vec &lik, int n );
	void condition( vec &cond );
};

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

template<class sq_T>
enorm<sq_T>::enorm( RV &rv, vec &mu0, sq_T &R0 ) {
	dim = rv.count();
	mu = mu0;
	R = R0;
};

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() {
	vec x( dim );
	RNG.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 ) {
	mat X( dim,N );
	vec x( dim );
	vec pom;
	int i;
	for ( i=0;i<N;i++ ) {
		RNG.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 ) {
	//
};


template<class sq_T>
enorm<sq_T>::enorm() {};

template<class sq_T>
mlnorm<sq_T>::mlnorm( RV &rv,RV &rvc, mat &A, sq_T &R ) {
	int dim = rv.count();
	vec mu( dim );

	epdf = enorm<sq_T>( rv,mu,R );
}

template<class sq_T>
vec mlnorm<sq_T>::samplecond( 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( 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<dim; 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( vec &cond ) {
	epdf.mu = A*cond;
//R is already assigned;
}

#endif //EF_H