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

//this comes first because it is used inside emix!

/*! \brief Class representing ratio of two densities
which arise e.g. by applying the Bayes rule.
It represents density in the form:
\f[
f(rv|rvc) = \frac{f(rv,rvc)}{f(rvc)}
\f]
where \f$ f(rvc) = \int f(rv,rvc) d\ rv \f$.

In particular this type of arise by conditioning of a mixture model.

At present the only supported operation is evalcond().
 */
class mratio: public mpdf {
	protected:
	//! Nominator in the form of mpdf
		const epdf* nom;
	//!Denominator in the form of epdf
		epdf* den;
	//!flag for destructor 
		bool destroynom;
	public:
	//!Default constructor. By default, the given epdf is not copied! 
	//! It is assumed that this function will be used only temporarily.
		mratio(const epdf* nom0, const RV &rv, bool copy=false):mpdf(rv,RV()){
			if (copy){
//			nom = nom0->_copy_();
				destroynom=true;
			} else {
				nom = nom0;
				destroynom = false;
			}
			rvc = nom->_rv().subt(rv);
			it_assert_debug(rvc.length()>0,"Makes no sense to use this object!");
			den = nom->marginal(rvc);
		};
		double evalcond(const vec &val, const vec &cond) {
			return exp(nom->evalpdflog(concat(val,cond)) - den->evalpdflog(cond));
		}
	//! Object takes ownership of nom and will destroy it
		void ownnom(){destroynom=true;}
	//! Default destructor
		~mratio(){delete den; if (destroynom) {delete nom;}}
};

/*!
* \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;
	//!Flag if owning Coms
	bool destroyComs;
public:
	//!Default constructor
	emix ( const RV &rv ) : epdf ( rv ) {};
	//! Set weights \c w and components \c Coms
	//!By default Coms are copied inside. \param copy can be set to false if Coms live externally. Use method ownComs() if Coms should be destroyed by the destructor.
	void set_parameters ( const vec &w, const Array<epdf*> &Coms, bool copy=true );

	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 ) * exp(Coms ( i )->evalpdflog ( val ));}
		return log ( sum );
	};
	vec evalpdflog_m ( const mat &Val ) const {
		vec x=ones(Val.cols());
		vec y(Val.cols());
		for (int i = 0; i < w.length(); i++ ) {
			y = w ( i )*exp(Coms ( i )->evalpdflog_m ( Val ) );
			elem_mult_inplace(y,x); //result is in x
		}
		return log(x);
	};
	mat evalpdflog_M ( const mat &Val ) const {
		mat X(w.length(), Val.cols());
		for (int i = 0; i < w.length(); i++ ) {
			X.set_row(i, w ( i )*exp(Coms ( i )->evalpdflog_m ( Val ) ));
		}
		return X;
	};

	emix* marginal ( const RV &rv ) const;
	mratio* condition ( const RV &rv ) const; //why not mratio!!

//Access methods
	//! returns a pointer to the internal mean value. Use with Care!
	vec& _w() {return w;}
	virtual ~emix() {if ( destroyComs ) {for ( int i=0;i<Coms.length();i++ ) {delete Coms ( i );}}}
	//! Auxiliary function for taking ownership of the Coms()
	void ownComs() {destroyComs=true;}
};

/*! \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 compositepdf, public mpdf {
protected:
	//! pointers to epdfs - shortcut to mpdfs()._epdf()
	Array<epdf*> epdfs;
	//! Indeces of rvc in common rvc
	Array<ivec> irvcs_rvc;
	//! Indeces of common rvc in rvcs
	Array<ivec> irvc_rvcs;
public:
	/*!\brief Constructor from list of mFacs,
	*/
	mprod ( Array<mpdf*> mFacs ) : compositepdf ( mFacs ), mpdf ( getrv ( true ),RV() ), epdfs ( n ), irvcs_rvc ( n ), irvc_rvcs ( n ) {
		setrvc ( rv,rvc );
		setindices ( rv );
		for ( int i = 0;i < n;i++ ) {
			// establish link between rvc and rvcs
			rvc.dataind ( mpdfs ( i )->_rvc(), irvc_rvcs ( i ), irvcs_rvc ( i ) );
		}

		for ( int i=0;i<n;i++ ) {
			epdfs ( i ) =& ( mpdfs ( i )->_epdf() );
		}
	};

	double evalcond ( const vec &val, const vec &cond ) {
		int i;
		double res = 0.0;
		vec condi;
		for ( i = n - 1;i >= 0;i-- ) {
			if ( ( irvcs_rvc ( i ).length() > 0 ) || ( irvcs_rv ( i ).length() > 0 ) ) {
				condi = zeros ( irvcs_rvc ( i ).length() +irvcs_rv ( i ).length() );
				//copy part of the rvc into condi
				set_subvector ( condi, irvcs_rvc ( i ), get_vec ( cond,irvc_rvcs ( i ) ) );
				//copy part of the rv into condi
				set_subvector ( condi, irvcs_rv ( i ), get_vec ( val,irv_rvcs ( i ) ) );
				mpdfs ( i )->condition ( condi );
			}
			// add logarithms
			res += epdfs ( i )->evalpdflog ( val ( rvsinrv ( i ) ) );
		}
		return exp ( res );
	}
	vec samplecond ( const vec &cond, double &ll ) {
		vec smp=zeros ( rv.count() );
		vec condi;
		vec smpi;
		ll = 0;
		for ( int i = ( n - 1 );i >= 0;i-- ) {
			if ( ( irvcs_rvc ( i ).length() > 0 ) || ( irvcs_rv ( i ).length() > 0 ) ) {
				condi = zeros ( irvcs_rv ( i ).length() + irvcs_rvc ( i ).length() );
				// copy data in condition
				set_subvector ( condi,irvcs_rvc ( i ), get_vec ( cond, irvc_rvcs ( i ) ) );
				// copy data in already generated sample
				set_subvector ( condi,irvcs_rv ( i ), get_vec ( smp,irv_rvcs ( i ) ) );

				mpdfs ( i )->condition ( condi );
			}
			smpi = epdfs ( i )->sample();
			// copy contribution of this pdf into smp
			set_subvector ( smp,rvsinrv ( i ), smpi );
			// add ith likelihood contribution
			ll+=epdfs ( i )->evalpdflog ( smpi );
		}
		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;
	}

	~mprod() {};
};

//! Product of independent epdfs. For dependent pdfs, use mprod.
class eprod: public epdf {
protected:
	//! Components (epdfs)
	Array<const epdf*> epdfs;
	//! Array of indeces
	Array<ivec> rvinds;
public:
	eprod ( const Array<const epdf*> epdfs0 ) : epdf ( RV() ),epdfs ( epdfs0 ),rvinds ( epdfs.length() ) {
		bool independent=true;
		for ( int i=0;i<epdfs.length();i++ ) {
			independent=rv.add ( epdfs ( i )->_rv() );
			it_assert_debug ( independent==true, "eprod:: given components are not independent ." );
			rvinds ( i ) = ( epdfs ( i )->_rv() ).dataind ( rv );
		}
	}

	vec mean() const {
		vec tmp ( rv.count() );
		for ( int i=0;i<epdfs.length();i++ ) {
			vec pom = epdfs ( i )->mean();
			set_subvector ( tmp,rvinds ( i ), pom );
		}
		return tmp;
	}
	vec sample() const {
		vec tmp ( rv.count() );
		for ( int i=0;i<epdfs.length();i++ ) {
			vec pom = epdfs ( i )->sample();
			set_subvector ( tmp,rvinds ( i ), pom );
		}
		return tmp;
	}
	double evalpdflog ( const vec &val ) const {
		double tmp=0;
		for ( int i=0;i<epdfs.length();i++ ) {
			tmp+=epdfs ( i )->evalpdflog ( val ( rvinds ( i ) ) );
		}
		return tmp;
	}
	//!access function
	const epdf* operator () ( int i ) const {it_assert_debug ( i<epdfs.length(),"wrong index" );return epdfs ( i );}
};


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