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

#define LOG2  0.69314718055995

#include "../shared_ptr.h"
#include "exp_family.h"

namespace bdm {

//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 evallogcond().
 */
class mratio: public pdf {
protected:
	//! Nominator in the form of pdf
	const epdf* nom;

	//!Denominator in the form of epdf
	shared_ptr<epdf> den;

	//!flag for destructor
	bool destroynom;
	//!datalink between conditional and nom
	datalink_m2e dl;
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 ) : pdf ( ), dl ( ) {
		// adjust rv and rvc

		set_rv ( rv ); 
		dim = rv._dsize();

		rvc = nom0->_rv().subt ( rv );
		dimc = rvc._dsize();

		//prepare data structures
		if ( copy ) {
			bdm_error ( "todo" );
			// destroynom = true;
		} else {
			nom = nom0;
			destroynom = false;
		}
		bdm_assert_debug ( rvc.length() > 0, "Makes no sense to use this object!" );

		// build denominator
		den = nom->marginal ( rvc );
		dl.set_connection ( rv, rvc, nom0->_rv() );
	};

	double evallogcond ( const vec &val, const vec &cond ) {
		double tmp;
		vec nom_val ( dimension() + dimc );
		dl.pushup_cond ( nom_val, val, cond );
		tmp = exp ( nom->evallog ( nom_val ) - den->evallog ( cond ) );
		return tmp;
	}

	//! Returns a sample from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv
	virtual vec samplecond ( const vec &cond ) NOT_IMPLEMENTED(0);

	//! Object takes ownership of nom and will destroy it
	void ownnom() {
		destroynom = true;
	}
	//! Default destructor
	~mratio() {
		if ( destroynom ) {
			delete nom;
		}
	}




private:
	// not implemented
	mratio ( const mratio & );
	mratio &operator= ( const mratio & );
};

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

public:
	//! Default constructor
	emix ( ) : epdf ( ) { }

	virtual void validate ();

	vec sample() const;

	vec mean() const;

	vec variance() const;

	double evallog ( const vec &val ) const;

	vec evallog_mat ( const mat &Val ) const;

	//! Auxiliary function that returns pdflog for each component
	mat evallog_coms ( const mat &Val ) const;

	shared_ptr<epdf> marginal ( const RV &rv ) const;
	//! Update already existing marginal density  \c target
	void marginal ( const RV &rv, emix &target ) const;
	shared_ptr<pdf> condition ( const RV &rv ) const;

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

	/*!
	  \brief returns a reference to the internal array of components. Use with Care! Set components \c Coms

	  Shared pointers in Coms are kept inside this instance and
	  shouldn't be modified after being passed to this method.
	*/
	Array<shared_ptr<epdf> >& _Coms ( ) {
		return Coms;
	}

	//! returns a reference to the internal components specified by index i. Use with Care!
	shared_ptr<epdf> _Coms ( int i ) {
		return Coms ( i );
	}

	void set_rv ( const RV &rv ) {
		epdf::set_rv ( rv );
		for ( int i = 0; i < Coms.length(); i++ ) {
			Coms ( i )->set_rv ( rv );
		}
	}

	//! Load from structure with elements:
	//!  \code
	//! { class='emix';
	//!   pdfs = (..., ...);     // list of pdfs in the mixture
	//!   weights = ( 0.5, 0.5 ); // weights of pdfs in the mixture
	//! }
	//! \endcode
	//!@}
	void from_setting ( const Setting &set );
};
SHAREDPTR ( emix );
UIREGISTER ( emix );

/*!
* \brief Mixture of egiws

*/
class egiwmix : public egiw {
protected:
	//! weights of the components
	vec w;
	//! Component (epdfs)
	Array<egiw*> Coms;
	//!Flag if owning Coms
	bool destroyComs;
public:
	//!Default constructor
	egiwmix ( ) : egiw ( ) {};

	//! Set weights \c w and components \c Coms
	//!By default Coms are copied inside. Parameter \c 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<egiw*> &Coms, bool copy = false );

	//!return expected value
	vec mean() const;

	//!return a sample from the density
	vec sample() const;

	//!return the expected variance
	vec variance() const;

	// TODO!!! Defined to follow ANSI and/or for future development
	void mean_mat ( mat &M, mat&R ) const {};
	double evallog_nn ( const vec &val ) const {
		return 0;
	};
	double lognc () const {
		return 0;
	}

	shared_ptr<epdf> marginal ( const RV &rv ) const;
	//! marginal density update
	void marginal ( const RV &rv, emix &target ) const;

//Access methods
	//! returns a pointer to the internal mean value. Use with Care!
	vec& _w() {
		return w;
	}
	virtual ~egiwmix() {
		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;
	}

	//!access function
	egiw* _Coms ( int i ) {
		return Coms ( i );
	}

	void set_rv ( const RV &rv ) {
		egiw::set_rv ( rv );
		for ( int i = 0; i < Coms.length(); i++ ) {
			Coms ( i )->set_rv ( rv );
		}
	}

	//! Approximation of a GiW mix by a single GiW pdf
	egiw* approx();
};

/*! \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 pdf {
private:
	Array<shared_ptr<pdf> > pdfs;

	//! Data link for each pdfs
	Array<shared_ptr<datalink_m2m> > dls;

public:
	//! \brief Default constructor
	mprod() { }

	/*!\brief Constructor from list of mFacs
	*/
	mprod ( const Array<shared_ptr<pdf> > &mFacs ) {
		set_elements ( mFacs );
	}
	//! Set internal \c pdfs from given values
	void set_elements ( const Array<shared_ptr<pdf> > &mFacs );

	double evallogcond ( const vec &val, const vec &cond );

	vec evallogcond_mat ( const mat &Dt, const vec &cond );

	vec evallogcond_mat ( const Array<vec> &Dt, const vec &cond );

	//TODO smarter...
	vec samplecond ( const vec &cond ) {
		//! Ugly hack to help to discover if mpfs are not in proper order. Correct solution = check that explicitely.
		vec smp = std::numeric_limits<double>::infinity() * ones ( dimension() );
		vec smpi;
		// Hard assumption here!!! We are going backwards, to assure that samples that are needed from smp are already generated!
		for ( int i = ( pdfs.length() - 1 ); i >= 0; i-- ) {
			// generate contribution of this pdf
			smpi = pdfs ( i )->samplecond ( dls ( i )->get_cond ( smp , cond ) );
			// copy contribution of this pdf into smp
			dls ( i )->pushup ( smp, smpi );
		}
		return smp;
	}

	//! Load from structure with elements:
	//!  \code
	//! { class='mprod';
	//!   pdfs = (..., ...);     // list of pdfs in the order of chain rule
	//! }
	//! \endcode
	//!@}
	void from_setting ( const Setting &set ) {
		Array<shared_ptr<pdf> > temp_array; 
		UI::get ( temp_array, set, "pdfs", UI::compulsory );
		set_elements ( temp_array );
	}
};
UIREGISTER ( mprod );
SHAREDPTR ( 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<datalink*> dls;
public:
	//! Default constructor
	eprod () : epdfs ( 0 ), dls ( 0 ) {};
	//! Set internal
	void set_parameters ( const Array<const epdf*> &epdfs0, bool named = true ) {
		epdfs = epdfs0;//.set_length ( epdfs0.length() );
		dls.set_length ( epdfs.length() );

		bool independent = true;
		if ( named ) {
			for ( int i = 0; i < epdfs.length(); i++ ) {
				independent = rv.add ( epdfs ( i )->_rv() );
				bdm_assert_debug ( independent, "eprod:: given components are not independent." );
			}
			dim = rv._dsize();
		} else {
			dim = 0;
			for ( int i = 0; i < epdfs.length(); i++ ) {
				dim += epdfs ( i )->dimension();
			}
		}
		//
		int cumdim = 0;
		int dimi = 0;
		int i;
		for ( i = 0; i < epdfs.length(); i++ ) {
			dls ( i ) = new datalink;
			if ( named ) {
				dls ( i )->set_connection ( epdfs ( i )->_rv() , rv );
			} else {
				dimi = epdfs ( i )->dimension();
				dls ( i )->set_connection ( dimi, dim, linspace ( cumdim, cumdim + dimi - 1 ) );
				cumdim += dimi;
			}
		}
	}

	vec mean() const;

	vec variance() const;

	vec sample() const;

	double evallog ( const vec &val ) const;

	//!access function
	const epdf* operator () ( int i ) const {
		bdm_assert_debug ( i < epdfs.length(), "wrong index" );
		return epdfs ( i );
	}

	//!Destructor
	~eprod() {
		for ( int i = 0; i < epdfs.length(); i++ ) {
			delete dls ( i );
		}
	}
};


/*! \brief Mixture of pdfs with constant weights, all pdfs are of equal RV and RVC

*/
class mmix : public pdf {
protected:
	//! Component (pdfs)
	Array<shared_ptr<pdf> > Coms;
	//!weights of the components
	vec w;
public:
	//!Default constructor
	mmix() : Coms ( 0 ) { }

	double evallogcond ( const vec &dt, const vec &cond ) {
		double ll = 0.0;
		for ( int i = 0; i < Coms.length(); i++ ) {
			ll += Coms ( i )->evallogcond ( dt, cond );
		}
		return ll;
	}

	vec samplecond ( const vec &cond );

	//! Load from structure with elements:
	//!  \code
	//! { class='mmix';
	//!   pdfs = (..., ...);     // list of pdfs in the mixture
	//!   weights = ( 0.5, 0.5 ); // weights of pdfs in the mixture
	//! }
	//! \endcode
	//!@}
	void from_setting ( const Setting &set );

	virtual void validate();
};
SHAREDPTR ( mmix );
UIREGISTER ( mmix );

}
#endif //MX_H