#include <stat/libEF.h>
#include <stat/emix.h>
using namespace itpp;

//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->evalpdflog ( rgr ) );
		}
	}
	return sumsum ( PPdf ) *xstep*ystep;
}

int main() {
	RNG_randomize();

	cout << "Testing enorm(2,2)"<<endl;
	// Setup model
	vec mu ( "1.1 -1" );
	ldmat R ( mat ( "1 -0.5; -0.5 2" ) );

	RV x ( "{x }" );
	RV y ( "{y }" );

	enorm<ldmat> E ( concat ( x,y ) );
	E.set_parameters ( mu,R );
	cout << "enorm mean value:" << E.mean() <<endl;
	cout << "enorm mean variance:" << R.to_mat() <<endl;
	cout << "enorm normalizing constant:" << E.lognc() <<endl;
	cout << " Numerically integrates to [1.0]: " << normcoef ( &E, vec ( "-5 5" ), vec ( "-5 5" ) ) << endl;

	int N=1000;
	vec ll ( N );
	mat Smp=E.sample ( 1000 );
	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;

	epdf* Mg = E.marginal ( y );
	mpdf* Cn = E.condition ( x );


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

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

	Smp=mEp.samplecond ( vec ( 0 ),ll,1000 );
	Emu = sum ( Smp,2 ) /N;
	Er = ( Smp*Smp.transpose() ) /N - outer_product ( Emu,Emu );

	cout << "composite mean: " <<Emu <<endl;
	cout << "composite variance: " <<Er <<endl;
	
	cout << endl << " test of pdflog at zero "<<endl;
	cout << "original: " << exp(E.evalpdflog(vec("0 0"))) <<endl;
	cout << "composite: " << mEp.evalcond(vec("0 0"),vec(0)) << endl;
}