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

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

void Test ( const egiw &E ) {
}

int main() {
	cout << "Testing eGiw(1,1)"<<endl;
	{
		// Setup model
		double mu=1.1;
		double s=0.1;

		// TEST 1x1 EGIW
		mat V ( 2,2 );
		V ( 0,0 ) = pow ( mu,2 ) +s;
		V ( 1,0 ) = mu;
		V ( 0,1 ) = V ( 1,0 );
		V ( 1,1 ) = 1.0;

		double nu=10;

		RV thr ( "{thr }", "2" );
		egiw E ( thr,nu*V,nu );
		cout << "egiw mean value:" << E.mean() <<endl;
		cout << "egiw normalizing constant:" << E.lognc() <<endl;

		int n=100;
		vec t_val ( 2 );

		mat pdf ( 2*n,n );
		vec Mu ( 2*n );
		vec Si ( n );

		for ( int i=0;i<2*n;i++ ) {
			Mu ( i ) = -2+i* ( 1.0/ ( n ) ) *3.0;
			t_val ( 0 ) = Mu ( i );
			for ( int j=0;j<n;j++ ) {
				Si ( j ) = ( j+1 ) * ( 1.0/n ) *2;
				t_val ( 1 ) = Si ( j );

				pdf ( i,j ) =E.evalpdflog ( t_val );
			}
		}

		mat Pdf=exp ( pdf );
		vec fm=sum ( Pdf,2 );
		vec fs=sum ( Pdf,1 );
		cout << "Numerical mean: " << vec_2 ( Mu*fm/sum ( fm ), Si*fs/sum ( fs ) ) <<endl;
		cout << "Numerical integral of pdf: "<<sumsum ( Pdf/n/n*3*2 ) <<endl;
	}
	cout << "Testing Egiw(1,2)"<<endl;
	{
		// Setup model
		double mu=1.1; //unit step parametr
		double b=3.0; // sequence of <1 -1 1 -1...>
		double s=0.1;


		// TEST 1x1 EGIW
		mat V ( 3,3 );
		V ( 0,0 ) = pow ( mu,2 ) +pow ( b ,2 )  +s;
		V ( 1,0 ) = mu;
		V ( 2,0 ) = b;

		V ( 0,1 ) = V ( 1,0 );
		V ( 1,1 ) = 1.0;
		V ( 2,1 ) = 0.0;

		V ( 0,2 ) = V ( 2,0 );
		V ( 1,2 ) = V ( 2,1 );
		V ( 2,2 ) = 1.0;


		double nu=20;

		RV thr ( "{thr }", "3" );
		egiw E ( thr,nu*V,nu );
		cout << "egiw mean value:" << E.mean() <<endl;
		cout << "egiw normalizing constant:" << E.lognc() <<endl;

		int n=100;
		vec t_val ( 3 );

		mat Tmp= zeros ( 2*n,n );

		double summ=0.0;
		for ( int k=0;k<n;k++ ) { // ALL b
			t_val ( 1 ) = 1 + k* ( 1.0/n ) * 4.0;
			for ( int i=0;i<2*n;i++ ) { //ALL mu
				t_val ( 0 ) = -2+i* ( 1.0/  n ) *3.0;
				for ( int j=0;j<n;j++ ) { // All sigma
					t_val ( 2 ) = ( j+1 ) * ( 1.0/n ) *2.0;

					Tmp ( i,j ) = E.evalpdflog ( t_val );
				}
			}
			summ += sumsum ( exp ( Tmp ) ) /n/n/n*3.0*2.0*4.0;
		}


		cout << "Numerical integral of pdf: "<<summ <<endl;
	}
	
	
}