/*!
  \file
  \brief Mergers for combination of pdfs
  \author Vaclav Smidl.

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

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

#ifndef MER_H
#define MER_H


#include "mixef.h"

namespace bdm{
using std::string;

/*!
@brief Function for general combination of pdfs

Mixtures of Gaussian densities are used internally. Switching to other densities should be trivial.
*/

class merger : public compositepdf, public epdf {
protected:
	//!Internal mixture of EF models
	MixEF Mix;
	//! Data link for each mpdf in mpdfs
	Array<datalink_m2e*> dls;
	//! Array of rvs that are not modelled by mpdfs at all (aux)
	Array<RV> rvzs;
	//! Data Links of rv0 mpdfs - these will be conditioned the (rv,rvc) of mpdfs
	Array<datalink_m2e*> zdls;

	//!Number of samples used in approximation
	int Ns;
	//!Number of components in a mixture
	int Nc;
	//!Prior on the log-normal merging model
	double beta;
	//! Projection to empirical density
	eEmp eSmp;

public:
//!Default constructor
	merger ( const Array<mpdf*> &S ) :
			compositepdf ( S ), epdf ( ),
			Mix ( Array<BMEF*> ( 0 ),vec ( 0 ) ), dls ( n ), rvzs ( n ), zdls ( n ), eSmp() {
		RV ztmp;
		// Extend rv by rvc!
		RV rvc; setrvc ( rv,rvc );
		rv.add ( rvc );
		for ( int i=0;i<n;i++ ) {
			//Establich connection between mpdfs and merger
			dls ( i ) = new datalink_m2e ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), rv );
			// find out what is missing in each mpdf
			ztmp= mpdfs ( i )->_rv();
			ztmp.add ( mpdfs ( i )->_rvc() );
			rvzs ( i ) =rv.subt ( ztmp );
			zdls ( i ) = new datalink_m2e ( rvzs ( i ), ztmp, rv ) ;
		};
		//Set Default values of parameters
		beta=2.0;
		Ns=100;
		Nc=10;
		Mix.set_method ( EM );
	}
//! Set internal parameters used in approximation
	void set_parameters ( double beta0, int Ns0, int Nc0 ) {beta=beta0;Ns=Ns0;Nc=Nc0;eSmp.set_n(Ns0,false);}
//!Initialize the proposal density. This function must be called before merge()!
	void init() { ////////////// NOT FINISHED
		Array<vec> Smps ( n );
		//Gibbs sampling
		for ( int i=0;i<n;i++ ) {Smps ( i ) =zeros ( 0 );}
	}
//!Create a mixture density using known proposal
	void merge ( const epdf* g0 );
//!Create a mixture density, make sure to call init() before the first call
	void merge () {merge ( & ( Mix._epdf() ) );};

//! Merge log-likelihood values
	vec lognorm_merge ( mat &lW );
//! sample from merged density
//! weight w is a
	vec sample ( ) const { return Mix._epdf().sample();}
	double evallog ( const vec &dt ) const {
		vec dtf=ones ( dt.length() +1 );
		dtf.set_subvector ( 0,dt );
		return Mix.logpred ( dtf );
	}
	vec mean() const {
		const Vec<double> &w = eSmp._w();
		const Array<vec> &S = eSmp._samples();
		vec tmp=zeros ( dim); 
		for ( int i=0; i<Ns; i++ ) {
			tmp+=w ( i ) *S ( i );
		}
		return tmp;
	}
	mat covariance() const {
		const vec &w = eSmp._w();
		const Array<vec> &S = eSmp._samples();
		
		vec mea = mean(); 
		
		cout << sum(w) << "," << w*w <<endl;
		
		mat Tmp=zeros(dim, dim);
		for ( int i=0; i<Ns; i++ ) {
			Tmp+=w ( i ) *outer_product(S ( i ), S(i));
		}
		return Tmp-outer_product(mea,mea);
	}
	vec variance() const {
		const vec &w = eSmp._w();
		const Array<vec> &S = eSmp._samples();
		
		vec tmp=zeros(dim);
		for ( int i=0; i<Ns; i++ ) {
			tmp+=w ( i ) *pow(S ( i ),2);
		}
		return tmp-pow(mean(),2);
	}
//! for future use
	virtual ~merger() {
		for ( int i=0; i<n; i++ ) {
			delete dls ( i );
			delete zdls ( i );
		}
	};

//! Access function
	MixEF& _Mix() {return Mix;}
//! Access function
	eEmp& _Smp() {return eSmp;}
};

}

#endif // MER_H
