/*!
  \file
  \brief Probability distributions for Mixtures of pdfs
  \author Vaclav Smidl.

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

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

#ifndef MX_H
#define MX_H

#include "libBM.h"
#include "libEF.h"
//#include <std>

using namespace itpp;

/*!
* \brief Mixture of epdfs

Density function:
\f[
f(x) = \sum_{i=1}^{n} w_{i} f_i(x), \quad \sum_{i=1}^n w_i = 1.
\f]
where \f$f_i(x)\f$ is any density on random variable \f$x\f$, called \a component,

*/
class emix : public epdf {
protected:
	//! weights of the components
	vec w;
	//! Component (epdfs)
	Array<epdf*> Coms;
public:
	//!Default constructor
	emix ( RV &rv ) : epdf ( rv ) {};
	//! Set weights \c w and components \c R
	void set_parameters ( const vec &w, const Array<epdf*> &Coms );

	vec sample() const;
	vec mean() const {
		int i; vec mu = zeros ( rv.count() );
		for ( i = 0;i < w.length();i++ ) {mu += w ( i ) * Coms ( i )->mean(); }
		return mu;
	}
	double evalpdflog ( const vec &val ) const {
		int i;
		double sum = 0.0;
		for ( i = 0;i < w.length();i++ ) {sum += w ( i ) * Coms ( i )->evalpdflog ( val );}
		return log ( sum );
	};

//Access methods
	//! returns a pointer to the internal mean value. Use with Care!
	vec& _w() {return w;}
};

/*! \brief Chain rule decomposition of epdf

Probability density in the form of Chain-rule decomposition:
\[
f(x_1,x_2,x_3) = f(x_1|x_2,x_3)f(x_2,x_3)f(x_3)
\]
Note that
*/
class mprod: public mpdf {
protected:
	//
	int n;
	// pointers to epdfs
	Array<epdf*> epdfs;
	Array<mpdf*> mpdfs;
	//
	//! Indeces of rvs in common rv
	Array<ivec> rvinds;
	//! Indeces of rvc in common rv
	Array<ivec> rvcinrv;
	//! Indeces of rvc in common rvc
	Array<ivec> rvcinds;
	//! Indicate independence of its factors
	bool independent;
//	//! Indicate internal creation of mpdfs which must be destroyed
//	bool intermpdfs;
public:
	//!Constructor from list of eFacs or list of mFacs
	mprod ( Array<mpdf*> mFacs ) : mpdf ( RV(), RV() ), n ( mFacs.length() ), epdfs ( n ), mpdfs ( mFacs ), rvinds ( n ), rvcinrv ( n ), rvcinds ( n ) {
		int i;
		// Create rv
		for ( i = 0;i < n;i++ ) {
			rv.add ( mpdfs ( i )->_rv() ); //add rv to common rvs.
			epdfs ( i ) = & ( mpdfs ( i )->_epdf() ); // add pointer to epdf
		};
		// Create rvc
		for ( i = 0;i < n;i++ ) {
			rvc.add ( mpdfs ( i )->_rv().subt ( rv ) ); //add rv to common rvs.
		};

		independent = true;
		//test rvc of mpdfs and fill rvinds
		for ( i = 0;i < n;i++ ) {
			// find ith rv in common rv
			rvinds ( i ) = mpdfs ( i )->_rv().dataind ( rv );
			// find ith rvc in common rv
			rvcinrv ( i ) = mpdfs ( i )->_rvc().dataind ( rv );
			// find ith rvc in common rv
			rvcinds ( i ) = mpdfs ( i )->_rvc().dataind ( rvc );
			//
			if ( rvcinds ( i ).length() >0 ) {independent = false;}
			if ( rvcinds ( i ).length() >0 ) {independent = false;}
		}
	};

	double evalpdflog ( const vec &val ) const {
		int i;
		double res = 0.0;
		for ( i = n - 1;i > 0;i++ ) {
			if ( rvcinds ( i ).length() > 0 )
				{mpdfs ( i )->condition ( val ( rvcinds ( i ) ) );}
			// add logarithms
			res += epdfs ( i )->evalpdflog ( val ( rvinds ( i ) ) );
		}
		return res;
	}
	vec samplecond ( const vec &cond, double &ll ) {
		vec smp=zeros ( rv.count() );
		vec condi;
		for ( int i = ( n - 1 );i >= 0;i-- ) {
			if ( rvcinds ( i ).length() > 0 ) {
				condi = zeros ( rvcinrv.length() + rvcinds.length() );
				// copy data in condition
				set_subvector ( condi,rvcinds ( i ), cond );
				// copy data in already generated sample
				set_subvector ( condi,rvcinrv ( i ), smp );

				mpdfs ( i )->condition ( condi );
			}
			// copy contribution of this pdf into smp
			set_subvector ( smp,rvinds ( i ), epdfs ( i )->sample() );
		}
		return smp;
	}
	mat samplecond ( const vec &cond, vec &ll, int N ) {
		mat Smp(rv.count(),N);
		for(int i=0;i<N;i++){Smp.set_col(i,samplecond(cond,ll(i)));}
		return Smp;
	}
// 	vec mean() const {
// 		vec tmp ( rv.count() );
// 		if ( independent ) {
// 			for ( int i=0;i<n;i++ ) {
// 				vec pom = epdfs ( i )->mean();
// 				set_subvector ( tmp,rvinds ( i ), pom );
// 			}
// 		}
// 		else {
// 			int N=50*rv.count();
// 			it_warning ( "eprod.mean() computed by sampling" );
// 			tmp = zeros ( rv.count() );
// 			for ( int i=0;i<N;i++ ) { tmp += sample();}
// 			tmp /=N;
// 		}
// 		return tmp;
// 	}

	~mprod() {};
};

/*! \brief Mixture of mpdfs with constant weights, all mpdfs are of equal type

*/
class mmix : public mpdf {
protected:
	//! Component (epdfs)
	Array<mpdf*> Coms;
	//!Internal epdf
	emix Epdf;
public:
	//!Default constructor
	mmix ( RV &rv, RV &rvc ) : mpdf ( rv, rvc ), Epdf ( rv ) {ep = &Epdf;};
	//! Set weights \c w and components \c R
	void set_parameters ( const vec &w, const Array<mpdf*> &Coms ) {
		Array<epdf*> Eps ( Coms.length() );

		for ( int i = 0;i < Coms.length();i++ ) {
			Eps ( i ) = & ( Coms ( i )->_epdf() );
		}
		Epdf.set_parameters ( w, Eps );
	};

	void condition ( const vec &cond ) {
		for ( int i = 0;i < Coms.length();i++ ) {Coms ( i )->condition ( cond );}
	};
};
#endif //MX_H