#include "epdf_harness.h"
#include "stat/exp_family.h"
#include "mat_checks.h"
#include "test_util.h"
#include "UnitTest++.h"

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());
    }
}

}

namespace bdm {

template<>
const ParticularUI<epdf_harness> &ParticularUI<epdf_harness>::factory(
    ParticularUI<epdf_harness>("epdf_harness"));

void epdf_harness::from_setting(const Setting &set) {
    hepdf = UI::build<epdf>(set, "epdf");
    UI::get(mean, set, "mean");
    UI::get(variance, set, "variance");

    if (set.exists("support")) {
        UI::get(support, set, "support");
    }

    if (set.exists("nsamples")) {
        UI::get(nsamples, set, "nsamples");
    }

    if (set.exists("integral")) {
        UI::get(integral, set, "integral");
    }

    if (set.exists("R")) {
        UI::get(R, set, "R");
    }
}

void epdf_harness::test()
{
    CHECK_CLOSE(mean, hepdf->mean(), epsilon);
    CHECK_CLOSE(variance, hepdf->variance(), epsilon);

    if (support.rows() == 2) {
	vec xb = support.get_row(0);
	vec yb = support.get_row(1);
	CHECK_CLOSE(mean, num_mean2(hepdf.get(), xb, yb, nsamples, nsamples), 0.1);
	CHECK_CLOSE(integral, normcoef(hepdf.get(), xb, yb, nsamples, nsamples), 0.1);
    }

    if (R.rows() > 0) {
	mat Smp = hepdf->sample_m(10000);
	int N = Smp.cols();
	vec Emu = Smp*ones(N) /N ;
	mat Er = (Smp*Smp.transpose())/N - outer_product(Emu,Emu);
	CHECK_CLOSE(mean, Emu, 0.1);
	CHECK_CLOSE(R, Er, 0.1);
    }
}

}