#include <itpp/itmex.h>
#include <estim/merger.h>
#include "config2mxstruct.h"
#include "mexlog.h"
#include "mexparse.h"

using namespace bdm;

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 Src;
	Src.addList(input[0],"Sources");
	
	UImxArray Sup (input[1]);
	UImxArray Mer (input[2]);
	// Sources
	Array<mpdf*> Sources;
	//abuse Mer to store sources

	Setting& _Sources=Src.lookup("Sources");
	int Slen=_Sources.getLength();
	Sources.set_size(Slen);
	for (int i=0; i<Slen; i++){
		mpdf* mtmp = UI::build<mpdf>(_Sources,i);
		if (mtmp) { 
			Sources(i)=mtmp;
		} else{ // source is not mpdf
			epdf* etmp = UI::build<epdf>(_Sources,i);
			if (etmp){
				Sources(i) = new mepdf(etmp, true);
			}
		}

	}

	merger_base* Merger=UI::build<merger_base>(Mer.getRoot());

	// Support
	Setting & _Supp=Sup.getRoot();
	
	Array<vec> bounds (0);
	UI::get (bounds, _Supp, "grid");
	ivec nbins(0);
	UI::get (nbins, _Supp, "nbins");
	
	epdf *g0=UI::build<epdf> (_Supp, "pdf");
	int npoints=100;
	_Supp.lookupValue("nsamples",npoints);
	
// COMPUTE RESULTS
	Merger->set_sources(Sources,true); // takes care of deletion of sources
	if (bounds.length() > 0 && nbins.length()>0) {
		Merger->set_support (bounds,nbins);
	} else {
		if (g0) {
			Merger->set_support (*g0,npoints);
			delete g0;
		}
	}
	Merger->merge();
	
	// Save results
	if (n_output>0){
		mxArray* 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);
		
		output[0] = tmp;
	}
}