#include "shared_ptr.h"
#include "stat/exp_family.h"
#include "stat/emix.h"

using namespace bdm;

//These lines are needed for use of cout and endl
using std::cout;
using std::endl;

double normcoef ( const epdf* ep,const vec &xb, const vec &yb, int Ngr=100 ) {
	mat PPdf ( Ngr+1,Ngr+1 );
	vec rgr ( 2 );

	int i=0,j=0;
	double xstep= ( xb ( 1 )-xb ( 0 ) ) /Ngr;
	double ystep= ( yb ( 1 )-yb ( 0 ) ) /Ngr;

	for ( double x=xb ( 0 );x<=xb ( 1 );x+= xstep,i++ ) {
		rgr ( 0 ) =x;j=0;
		for ( double y=yb ( 0 );y<=yb ( 1 );y+=ystep,j++ ) {
			rgr ( 1 ) =y;
			PPdf ( i,j ) =exp ( ep->evallog ( rgr ) );
		}
	}
	return sumsum ( PPdf ) *xstep*ystep;
}

int main() {
	RNG_randomize();

	RV x ( "{x }" );
	RV y ( "{y }" );
	RV xy=concat(x,y);
	
	enorm<ldmat> E1; E1.set_rv ( xy );
		
	E1.set_parameters( "1.00054 1.0455" , mat ( "0.740142 -0.259015; -0.259015 1.0302" ));
	enorm<ldmat> E2;E2.set_rv ( xy );
	E2.set_parameters( "-1.2 -0.1" , mat ( "1 0.4; 0.4 0.5" ));

	Array<epdf*> A1(1);
	A1(0) = &E1;
	
	emix M1; M1.set_rv(xy);
	M1.set_parameters(vec("1"), A1, false);
	cout << "======== test if ARX and emix with one ARX are the same ==="<<endl;
	
	epdf* Mm = M1.marginal(y);
	epdf* Am = E1.marginal(y);
	mpdf* Mc = M1.condition(y);
	mpdf* Ac = E1.condition(y);
	
	cout << *((mlnorm<ldmat>*)Ac) <<endl;
	
	cout << "Mix marg at 0: " << Mm->evallog(vec_1(0.0)) <<endl;
	cout << "ARX marg at 0: " << Am->evallog(vec_1(0.0)) <<endl;
	cout << "Mix cond at 0,0: " << Mc->evallogcond(vec_1(0.0),vec_1(0.0)) <<endl;
	cout << "ARX cond at 0,0: " << Ac->evallogcond(vec_1(0.0),vec_1(0.0)) <<endl;
	
	cout << "======== Mixture with two components ==="<<endl;
	Array<epdf*> A2(2);
	A2(0) = &E1;
	A2(1) = &E2;
	
	emix M2; M2.set_rv(xy);
	M2.set_parameters(vec("1"), A2, false);
	
	
	cout << "Mixture normalization: " << normcoef(&M2, vec("-3 3 "), vec("-3 3 ")) <<endl;
	
	int N=3;
	vec ll ( N );
	vec ll2 ( N );
	mat Smp=M2.sample_m ( N );
	ll = M2.evallog_m(Smp);
	
	vec Emu = sum ( Smp,2 ) /N;
	mat Er = ( Smp*Smp.transpose() ) /N - outer_product ( Emu,Emu );

	cout << "original empirical mean: " <<Emu <<endl;
	cout << "original empirical variance: " <<Er <<endl;

	shared_ptr<epdf> Mg(dynamic_cast<epdf *>(M2.marginal(y)));
	mpdf* Cn = (mpdf*)M2.condition ( x );

	it_file it("emix_test.it");
	it << Name("Smp") << Smp;

	cout<< "marginal mean: " << Mg->mean() <<endl;

	cout << "========== putting them back together ======= "<<endl;
	mepdf mMg ( Mg );
	Array<mpdf*> AA ( 2 );
	AA ( 0 ) = Cn;
	AA ( 1 ) = &mMg;
	mprod mEp ( AA );

	for (int j=0; j<N; j++){
		ll2(j)=mEp.evallogcond (Smp.get_col(j), vec ( 0 ));
	}
	it << Name("ll") << ll;
	it << Name("ll2") << ll2;
	
	
}