#define BDMLIB // not an ideal way to prevent double registration of UI factories...
#include "stat/exp_family.h"
#include "stat/emix.h"
#include "mat_checks.h"
#include "UnitTest++.h"

using namespace bdm;

const double epsilon = 0.00001;

namespace UnitTest
{

inline void CheckClose(TestResults &results, const itpp::vec &expected,
		       const itpp::vec &actual, double tolerance,
		       TestDetails const& details) {
    if (!AreClose(expected, actual, tolerance)) { 
        MemoryOutStream stream;
        stream << "Expected " << expected << " +/- " << tolerance << " but was " << actual;

        results.OnTestFailure(details, stream.GetText());
    }
}

inline void CheckClose(TestResults &results, const itpp::mat &expected,
		       const itpp::mat &actual, double tolerance,
		       TestDetails const& details) {
    if (!AreClose(expected, actual, tolerance)) { 
        MemoryOutStream stream;
        stream << "Expected " << expected << " +/- " << tolerance << " but was " << actual;

        results.OnTestFailure(details, stream.GetText());
    }
}

}

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

TEST(test_enorm) {
    RNG_randomize();

    // 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;
    E.set_rv(concat(x,y));
    E.set_parameters(mu, R);
    CHECK_EQUAL(mu, E.mean());
    CHECK_CLOSE(2.11768, E.lognc(), epsilon);
    CHECK_CLOSE(1.0, normcoef(&E, vec("-5 5"), vec("-5 5")), 0.01);

    int N = 1000;
    vec ll(N);
    mat Smp = E.sample(1000);
    vec Emu = sum(Smp, 2) / N;
    CHECK_CLOSE(mu, Emu, 0.3);

    mat Er = (Smp * Smp.T()) / N - outer_product(Emu, Emu);
    CHECK_CLOSE(R.to_mat(), Er, 0.3);

    epdf *Mg = E.marginal(y);
    CHECK_CLOSE(vec("-1"), Mg->mean(), epsilon);

    // putting them back together
    mpdf *Cn = E.condition(x);
    mepdf mMg(Mg);
    Array<mpdf *> A(2);
    A(0) = Cn;
    A(1) = &mMg;
    mprod mEp(A);
    Smp = mEp.samplecond(vec(0), 1000);
    Emu = sum(Smp, 2) / N;
    CHECK_CLOSE(mu, Emu, 0.3);

    Er = (Smp * Smp.T()) / N - outer_product(Emu, Emu);
    CHECK_CLOSE(R.to_mat(), Er, 0.3);

    // test of pdflog at zero
    vec zero(0);
    vec zero2("0 0");
    CHECK_CLOSE(E.evallog(zero2), mEp.evallogcond(zero2, zero), epsilon);
}

// from testEpdf
TEST(test_enorm_sum) {
    vec x = "-10:0.1:10";
    vec y = "-10:0.1:10";

    RV rv("{x2 }", "2");
    vec mu0 = "0.0 0.0";
    mat V0 = "5 -0.05; -0.05 5.20";
    fsqmat R(V0);

    enorm<fsqmat> eN;
    eN.set_rv(rv);
    eN.set_parameters(mu0, R);

    vec pom(2);
    double suma = 0.0;
    for (int i = 0; i < x.length(); i++) {
	for (int j=0; j<y.length(); j++) {
	    pom(0) = x(i);
	    pom(1) = y(j);
	    suma += exp(eN.evallog(pom));
	}
    }

    CHECK_CLOSE(100, suma, 0.1);
}