#include "mpdf_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<mpdf_harness> &ParticularUI<mpdf_harness>::factory(
    ParticularUI<mpdf_harness>("mpdf_harness"));

void mpdf_harness::from_setting(const Setting &set) {
    hmpdf = UI::build<mpdf>(set, "mpdf");
    UI::get(cond, set, "cond");
    UI::get(mean, set, "mean");

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

    UI::get(R, set, "R");
}

void mpdf_harness::test()
{
    mat smp = hmpdf->samplecond_m(cond, nsamples);
    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);
}

}