#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;
	shared_ptr<enorm<ldmat> > eN = new enorm<ldmat>();
	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;
	shared_ptr<egamma> eG = new egamma();
	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;
	shared_ptr<mgamma> mG = new mgamma();
	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;
	shared_ptr<emix> eMix = new emix();
	Array<shared_ptr<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 );

	// emix::set_parameters requires the first mpdf to be named
	mG->set_rv(x);
	mG->set_rvc(y);
	mComs ( 0 ) = mG;

	eN->set_mu ( vec_2 ( 0.0, 0.0 ) );
	shared_ptr<mepdf> mEnorm = new mepdf ( 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<shared_ptr<mpdf> > A ( 2 );
	A ( 0 ) = new mepdf ( eN );
	A ( 1 ) = new mepdf ( eG );

	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;

}