#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;

void disp ( const vec &tmu, const mat &tR, const mat &Smp ) {
	int N = Smp.cols();
	vec Emu = Smp * ones ( N ) / N ;
	mat Er = ( Smp * Smp.transpose() ) / N - outer_product ( Emu, Emu );
	cout << "True mu:" << tmu << endl;
	cout << "Emp  mu:" << Emu << endl;

	cout << "True R:" << tR << endl;
	cout << "Emp  R:" << Er << endl;
}

int main() {

	RNG_randomize();

	RV x ( "{x }", "2" );
	RV y ( "{y }", "2" );
	int N = 10000; //number of samples
	vec mu0 = "1.5 1.7";
	mat V0 ( "1.2 0.3; 0.3 5" );
	ldmat R = ldmat ( V0 );

	cout << "====== ENorm ====== " << endl;
	enorm<ldmat> eN;
	eN.set_parameters ( mu0, R );
	mat Smp = eN.sample_m ( N );

	disp ( mu0, R.to_mat(), Smp );

	cout << "====== MlNorm ====== " << endl;
	mat I = eye ( 2 );
	mlnorm<ldmat> ML;
	ML.set_parameters ( I, zeros ( 2 ), R );
	Smp = ML.samplecond_m ( mu0, N );

	disp ( mu0, R.to_mat(), Smp );

	cout << "====== EGamma ====== " << endl;
	vec a = "100000,10000";
	vec b = a / 10.0;
	egamma eG;
	eG.set_parameters ( a, b );

	cout << eG.evallog ( a ) << endl;
	Smp = eG.sample_m ( N );

	vec g_mu = elem_div ( a, b );
	vec g_var = elem_div ( a, pow ( b, 2.0 ) );
	disp ( g_mu, diag ( g_var ), Smp );

	cout << "====== MGamma ====== " << endl;
	mgamma mG;
	double k = 10.0;
	mG.set_parameters ( k, mu0 );

	Smp = mG.samplecond_m ( mu0, N );
	disp ( mu0, pow ( mu0, 2.0 ) / k, Smp );

	cout << "======= EMix ======== " << endl;
	emix eMix;
	Array<epdf*> Coms ( 2 );
	Coms ( 0 ) = &eG;
	Coms ( 1 ) = &eN;

	eMix.set_parameters ( vec_2 ( 0.5, 0.5 ), Coms );
	vec smp = eMix.sample();
	Smp = eMix.sample_m ( N );
	disp ( eMix.mean(), zeros ( 2 ), Smp );

	cout << "======= MEpdf ======== " << endl;
	mepdf meMix ( &eMix );

	Smp = meMix.samplecond_m ( mu0, N );
	disp ( eMix.mean(), zeros ( 2 ), Smp );

	cout << "======= MMix ======== " << endl;
	mmix mMix;
	Array<shared_ptr<mpdf> > mComs ( 2 );
	mComs ( 0 ) = &mG;
	eN.set_mu ( vec_2 ( 0.0, 0.0 ) );
	mepdf mEnorm ( &eN );
	mComs ( 1 ) = &mEnorm;
	mMix.set_parameters ( vec_2 ( 0.5, 0.5 ), mComs );

	Smp = mMix.samplecond_m ( mu0, N );
	disp ( 0.5*eN.mean()+0.4*eG.mean(), zeros ( 2 ), Smp );

	cout << "======= EProd ======== " << endl;
	// we have to change eG.rv to y
	eN.set_rv ( x );
	eG.set_rv ( y );
	//create array
	Array<mpdf*> A ( 2 );
	mepdf meN ( &eN );
	mepdf meG ( &eG );
	A ( 0 ) = &meN;
	A ( 1 ) = &meG;

	mprod eP ( A );
	mat epV = zeros ( 4, 4 );
	epV.set_submatrix ( 0, 0, V0 );
	epV.set_submatrix ( 2, 2, diag ( g_var ) );

	vec v0 = vec ( 0 );
	Smp = eP.samplecond ( v0, N );
	disp ( concat ( eN.mean(), eG.mean() ), epV, Smp );

	cout << "======= eWishart ======== " << endl;
	mat wM = "1.0 0.9; 0.9 1.0";
	eWishartCh eW;
	eW.set_parameters ( wM / 100, 100 );
	mat mea = zeros ( 2, 2 );
	mat Ch ( 2, 2 );
	for ( int i = 0; i < 100; i++ ) {
		Ch = eW.sample_mat();
		mea += Ch.T() * Ch;
	}
	cout << mea / 100 << endl;

	cout << "======= rwiWishart ======== " << endl;
	rwiWishartCh rwW;
	rwW.set_parameters ( 2, 0.1, "1 1", 0.9 );
	mea = zeros ( 2, 2 );
	mat wMch = chol ( wM );
	for ( int i = 0; i < 100; i++ ) {
		vec tmp = rwW.samplecond ( vec ( wMch._data(), 4 ) );
		copy_vector ( 4, tmp._data(), Ch._data() );
		mea += Ch.T() * Ch;
	}
	cout << mea / 100 << endl;
	//Exit program:
	return 0;

}