#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(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<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(mMix._epdf().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;

}