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


#include "../estim/mixtures.h"

namespace bdm
{
using std::string;

//!Merging methods
enum MERGER_METHOD {ARITHMETIC = 1, GEOMETRIC = 2, LOGNORMAL = 3};

/*!
@brief Base class for general combination of pdfs on discrete support

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

The merged pdfs are expected to be of the form:
 \f[ f(x_i|y_i),  i=1..n \f]
where the resulting merger is a density on \f$ \cup [x_i,y_i] \f$ .
Note that all variables will be joined.

As a result of this feature, each source must be extended to common support
\f[ f(z_i|y_i,x_i) f(x_i|y_i) f(y_i)  i=1..n \f]
where \f$ z_i \f$ accumulate variables that were not in the original source.
These extensions are calculated on-the-fly.

However, these operations can not be performed in general. Hence, this class merges only sources on common support, \f$ y_i={}, z_i={}, \forall i \f$.
For merging of more general cases, use offsprings merger_mix and merger_grid.
*/

class merger_base : public compositepdf, public epdf
{
	protected:
		//! Data link for each mpdf in mpdfs
		Array<datalink_m2e*> dls;
		//! Array of rvs that are not modelled by mpdfs at all, \f$ z_i \f$
		Array<RV> rvzs;
		//! Data Links for extension \f$ f(z_i|x_i,y_i) \f$
		Array<datalink_m2e*> zdls;
		//! number of support points
		int Npoints;
		//! number of sources
		int Nsources;

		//! switch of the methoh used for merging
		MERGER_METHOD METHOD;
		//!Prior on the log-normal merging model
		double beta;

		//! Projection to empirical density (could also be piece-wise linear)
		eEmp eSmp;

		//! debug or not debug
		bool DBG;

		//! debugging file
		it_file* dbg_file;
	public:
		//! \name Constructors
		//! @{

		//!Empty constructor
		merger_base () : compositepdf() {DBG=false;dbg_file=NULL;};
		//!Constructor from sources
		merger_base (const Array<mpdf*> &S) {set_sources (S);};
		//! Function setting the main internal structures 
		void set_sources (const Array<mpdf*> &Sources) {
			compositepdf::set_elements (Sources);
			//set sizes
			dls.set_size (Sources.length());
			rvzs.set_size (Sources.length());
			zdls.set_size (Sources.length());

			rv = getrv (/* checkoverlap = */ false);
			RV rvc; setrvc (rv, rvc);  // Extend rv by rvc!
			// join rv and rvc - see descriprion
			rv.add (rvc);
			// get dimension
			dim = rv._dsize();

			// create links between sources and common rv
			RV xytmp;
			for (int i = 0;i < mpdfs.length();i++) {
				//Establich connection between mpdfs and merger
				dls (i) = new datalink_m2e;
				dls (i)->set_connection (mpdfs (i)->_rv(), mpdfs (i)->_rvc(), rv);

				// find out what is missing in each mpdf
				xytmp = mpdfs (i)->_rv();
				xytmp.add (mpdfs (i)->_rvc());
				// z_i = common_rv-xy
				rvzs (i) = rv.subt (xytmp);
				//establish connection between extension (z_i|x,y)s and common rv
				zdls (i) = new datalink_m2e; zdls (i)->set_connection (rvzs (i), xytmp, rv) ;
			};
		}
		//! set debug file
		void set_debug_file (const string fname) { 
			if (DBG) delete dbg_file; 
			dbg_file = new it_file (fname); 
			if (dbg_file) DBG = true;
		}
		//! Set internal parameters used in approximation
		void set_method (MERGER_METHOD MTH, double beta0=0.0) {
			METHOD = MTH; 
			beta = beta0;
		}
		//! Set support points from a pdf by drawing N samples
		void set_support(const epdf &overall, int N){
			it_assert_debug ( rv.equal ( overall._rv() ),"Incompatible parameter overall!" );
			eSmp.set_statistics(&overall,N);
			Npoints=N;
		}
		
		//! Destructor
		virtual ~merger_base() {
			for (int i = 0; i < Nsources; i++) {
				delete dls (i);
				delete zdls (i);
			}
			if (DBG) delete dbg_file;
		};
		//!@}
		
		//! \name Mathematical operations
		//!@{
		
		//!Merge given sources in given points
		void merge () {
			if (eSmp._w().length() ==0) {it_error("Empty support points use set_support" );}
			//check if sources overlap:
			bool OK = true;
			for (int i = 0;i < mpdfs.length(); i++) {
				OK &= (rvzs (i)._dsize() == 0); // z_i is empty
				OK &= (mpdfs (i)->_rvc()._dsize() == 0); // y_i is empty
			}

			if (OK) {
				mat lW = zeros (mpdfs.length(), eSmp._w().length());

				vec emptyvec (0);
				for (int i = 0; i < mpdfs.length(); i++) {
					for (int j = 0; j < eSmp._w().length(); j++) {
						lW (i, j) = mpdfs (i)->evallogcond (eSmp._samples() (j), emptyvec);
					}
				}

				vec wtmp = exp (merge_points (lW));
				//renormalize
				eSmp._w() = wtmp / sum (wtmp);
			} else {
				it_error("Sources are not compatible - use merger_mix");
			}
			;
		};


		//! Merge log-likelihood values in points using method specified by parameter METHOD
		vec merge_points (mat &lW);
		
		
		//! sample from merged density
//! weight w is a
		vec mean() const {
			const Vec<double> &w = eSmp._w();
			const Array<vec> &S = eSmp._samples();
			vec tmp = zeros (dim);
			for (int i = 0; i < Npoints; 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 < Npoints; 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 < Nsources; i++) {
				tmp += w (i) * pow (S (i), 2);
			}
			return tmp -pow (mean(), 2);
		}
		//!@}

		//! \name Access to attributes
		//! @{

		//! Access function
		eEmp& _Smp() {return eSmp;}
		
		//!@}
};

class merger_mix : public merger_base
{
	protected:
		//!Internal mixture of EF models
		MixEF Mix;
		//!Number of components in a mixture
		int Ncoms;
		//! coefficient of resampling
		double effss_coef;

	public:
		//!\name Constructors
		//!@{
		merger_mix () {};
		merger_mix (const Array<mpdf*> &S) {set_sources(S);};
		//! Set sources and prepare all internal structures
		void set_sources (const Array<mpdf*> &S) {
			merger_base::set_sources(S);
			Nsources = S.length();
		}
		//! Set internal parameters used in approximation
		void set_parameters (int Ncoms0=10, double effss_coef0 = 0.5) {
			Ncoms = Ncoms0;
			effss_coef = effss_coef0;
		}
		//!@}
		
		//! \name Mathematical operations
		//!@{
		
		//!Merge values using mixture approximation
		void merge ();

		//! sample from the approximating mixture
		vec sample () const { return Mix.posterior().sample();}
		//! loglikelihood computed on mixture models
		double evallog (const vec &dt) const {
			vec dtf = ones (dt.length() + 1);
			dtf.set_subvector (0, dt);
			return Mix.logpred (dtf);
		}
		//!@}

		//!\name Access functions
		//!@{
//! Access function
		MixEF& _Mix() {return Mix;}
		//! Access function
		emix* proposal() {emix* tmp = Mix.epredictor(); tmp->set_rv (rv); return tmp;}
	//! @}
		
};

}

#endif // MER_H
