#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." );
	RV::clear_all();
	// 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<shared_ptr<mpdf> >  Sources;
	UI::get(Sources, Cfg, "Sources", UI::compulsory);
	shared_ptr<merger_base> Merger = UI::build<merger_base> ( Cfg, "Merger" );

	// Support
	try{
	shared_ptr<rectangular_support> RecSup = UI::build<rectangular_support> ( Cfg, "Support" ,UI::optional);
	if(RecSup){
		Merger->set_support ( *RecSup );
	} 
	} catch (UIException &e) {
		shared_ptr<discrete_support> DisSup = UI::build<discrete_support> ( Cfg, "Support", UI::optional );
		if (DisSup){
			Merger->set_support ( *DisSup );
		} else{bdm_error("Support not defined");}
	}
// COMPUTE RESULTS
	Merger->set_sources ( 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 );
	}

	merger_mix* MerMix=dynamic_cast<merger_mix*> ( Merger.get() );
	vec mix_val;

	if ( MerMix ) {
		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 );

		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
		if ( mix_val.length() >0 ) {
			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
}