#include <stat/merger.h>

#ifdef MEX
#include <itpp/itmex.h>
#include <mex/mex_logger.h>
#include <mex/mex_parser.h>
#include <mex/config2mxstruct.h>
#endif

/*! \file 
\brief Merging static pdfs
Merger is a program/mex-function for static merging of given pdfs on given support.

In command line, it expected config file with the following structure:
\code
Sources = (... any epdf or mpdf  ...);
Support = {grid=([low1d,high1d],[low2d,high2d],...);
		   nbins=[bins1d,bins2d,... ]}; 
// OR
Support = {pdf={class='epdf',...};
		   nsamples=2}; 
Merger = {class="merger_base",...}
\endcode

In mex-file, the above structures (Sources,Support,Merger) are input arguments. Type >>merger for help.
*/

using namespace bdm;

#ifdef MEX 
void mexFunction(int n_output, mxArray *output[], int n_input, const mxArray *input[])
{
    // Check the number of inputs and output arguments
	if(n_input!=3) mexErrMsgTxt("Usage:\n"  
		"result=merger(sources, support, merger)\n"
		"  sources= { struct('class','epdf'),... };  % cell of pdfs (epdfs or mpdfs) to be merged,\n"
		"  support= struct(\n"
		"           grid    = {[dim1_start,dim1_end], [dim2_start, dim2_end]...}  %support boundary \n"   
		"           nbins   = [bins_in_dim1, bins_in_dim2,...]                    %fixed \n"   
		"         === OR ==\n"
		"           pdf     = struct('class','epdf'); % pdf to draw samples from\n"
		"           nsamples= 100;                    % number of samples\n"
		"           );\n"
		"        If all elements are present,  (grid,nbins) is used;\n"
		"  merger = struct('class','merger_*');       % object to be used for merging,\n\n"
		"see documentation of classes epdf, mpdf, merger_base and their offsprings in BDM.");

	// LOAD CONFIG
	UImxArray Cfg;
	Cfg.addList(input[0],"Sources");
	Cfg.addGroup(input[1],"Support");
	Cfg.addGroup(input[2],"Merger");

	//DBG
	Cfg.writeFile("merger.cfg");
#else
	int main()
	{
		UIFile Cfg("merger.cfg");
#endif	
	// Sources
	Array<mpdf*> Sources;
	//abuse Mer to store sources
	Setting& _Sources=Cfg.lookup("Sources");
	int Slen=_Sources.getLength();
	Sources.set_size(Slen);
	for (int i=0; i<Slen; i++){
		try{
			mpdf* mtmp = UI::build<mpdf>(_Sources,i);
			Sources(i)=mtmp;
		}
		catch (UIException){
			// it is not mpdf - see if it is epdf
			try {
				epdf* etmp = UI::build<epdf>(_Sources,i);
				if (etmp){
					Sources(i) = new mepdf(etmp, true);
				}				
			}
			catch (UIException e) 
			{
				it_error("No mpdfs or epdfs found! " + string(e.what()));
			}
			catch (std::exception e) {
				it_error("Error in UI at "+_Sources[i].getPath());
			}
		}
		catch (std::exception e) {
			it_error("Error in UI at "+_Sources[i].getPath());
		}
	}

	merger_base* Merger=UI::build<merger_base>(Cfg,"Merger");

	// Support
	Setting & _Supp=Cfg.lookup("Support");
	
	if (_Supp.exists("grid") &&  _Supp.exists("nbins")) {
	Array<vec> bounds (0);
	UI::get (bounds, _Supp, "grid");
	ivec nbins(0);
	UI::get (nbins, _Supp, "nbins");
	Merger->set_support (bounds,nbins);
	
	}else {
		if (_Supp.exists("pdf") &&  _Supp.exists("nsamples")){
			epdf *g0=UI::build<epdf> (_Supp, "pdf");
			int npoints=100;
			_Supp.lookupValue("nsamples",npoints);
			Merger->set_support (*g0,npoints);
			delete g0;	
		}
		else it_error("Use either [grid,nbins] or [pdf,nsamples].");
	}
// COMPUTE RESULTS
	Merger->set_sources(Sources,true); // takes care of deletion of sources
	Merger->merge();
	
// save results
	Array<vec> source_vals(Sources.length());
		for (int i=0;i<Sources.length();i++){
			datalink_m2e dl;
			dl.set_connection(Sources(i)->_rv(), Sources(i)->_rvc(), Merger->_rv());
			
			vec ll(Merger->_Smp()._samples().length());
			for(int j=0; j<Merger->_Smp()._samples().length(); j++){
				vec dt=dl.pushdown(Merger->_Smp()._samples()(j));
				vec dtc=dl.get_cond(Merger->_Smp()._samples()(j));
				ll(j)=Sources(i)->evallogcond(dt,dtc);
			}

			vec sll = exp(ll);
			
			source_vals(i)=sll/sum(sll);
		}		

			vec ll(Merger->_Smp()._samples().length());
			for(int j=0; j<Merger->_Smp()._samples().length(); j++){
				ll(j)=Merger->evallog(Merger->_Smp()._samples()(j));
			}

			vec sll = exp(ll);
			
			vec mix_val=sll/sum(sll);
		
	
#ifdef MEX
	mxArray* tmp ;
	// Save results
	if (n_output>0){
		tmp = mxCreateStructMatrix(1,1,0,NULL);
		//support
		Array<vec> &samples=Merger->_Smp()._samples();
		if (samples.size()>0){
			mxArray* fld=mxCreateDoubleMatrix(samples(0).length(), samples.size(), mxREAL);
		Arrayvec2mxArray(samples,fld);
		mxReplaceFieldNM(tmp, "support", fld);
		}

		//weights
		vec &w = Merger->_Smp()._w();
		mxArray* fldw=mxCreateDoubleMatrix(1, w.length(), mxREAL);
		vec2mxArray(w,fldw);
		mxReplaceFieldNM(tmp, "weights", fldw);
		
		//mixture values
		mxArray* fldm=mxCreateDoubleMatrix(1, w.length(), mxREAL);
		vec2mxArray(mix_val,fldm);
		mxReplaceFieldNM(tmp, "mix", fldm);

		// sources
		char srcstr[20];
		for (int i=0;i<Sources.length();i++){
			sprintf(srcstr,"source%d",i+1);
			mxArray* fldw=mxCreateDoubleMatrix(1, source_vals(i).length(), mxREAL);
			vec2mxArray(source_vals(i),fldw);
			mxReplaceFieldNM(tmp, srcstr, fldw);
		}		

		output[0] = tmp;
	}
#endif
}