#include "epdf_harness.h"
#include "base/bdmbase.h"
#include "base/user_info.h"
#include "stat/exp_family.h"
#include "stat/emix.h"
#include "mat_checks.h"
#include "test_util.h"
#include "UnitTest++.h"
#include <memory>

namespace bdm {

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

void epdf_harness::test_config(const char *config_file_name) {
    RV::clear_all();

    UIFile in(config_file_name);
    Array<epdf_harness *> input;
    UI::get(input, in, "data");
    int sz = input.size();
    CHECK(sz > 0);
    for (int i = 0; i < sz; ++i) {
        input(i)->test(config_file_name, i);
    }
}

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


	UI::get(support, set, "support");
	UI::get(nbins, set, "nbins");
    UI::get(nsamples, set, "nsamples");
    UI::get(R, set, "R");

	RV* rv = UI::build<RV>(set, "marginal_rv");
    if (rv) 
        mrv = shared_ptr<RV>(rv);    

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

void epdf_harness::test(const char *config_name, int idx)
{
    CurrentContext cc(config_name, idx);

    CHECK_CLOSE_EX(mean, hepdf->mean(), tolerance);
    CHECK_CLOSE_EX(variance, hepdf->variance(), tolerance);

    if (support.rows() == 2) {
	vec xb = support.get_row(0);
	vec yb = support.get_row(1);

	int old_size = nbins.size();
	if (old_size < 2) {
	    vec new_nbins("100 100");
	    for (int i = 0; i < old_size; ++i) {
	        new_nbins(i) = nbins(i);
	    }

	    nbins = new_nbins;
	}

	CHECK_CLOSE_EX(mean, num_mean2(hepdf.get(), xb, yb, nbins(0), nbins(1)), tolerance);
	CHECK_CLOSE_EX(1.0, normcoef(hepdf.get(), xb, yb, nbins(0), nbins(1)), tolerance);
    }

    if (R.rows() > 0) {
	mat smp = hepdf->sample_m(nsamples);
	vec emu = smp * ones(nsamples) / nsamples;
	mat er = (smp * smp.T()) / nsamples - outer_product(emu, emu);

	// simplify overloading for Visual Studio
	vec delta = sqrt(variance) / sqrt(static_cast<double>(nsamples));
	CHECK_CLOSE_EX(mean, emu, delta);

	CHECK_CLOSE_EX(R, er, tolerance);
    }

    if (mrv.get()) {
        RV crv = hepdf->_rv().subt(*mrv);
	shared_ptr<epdf> m = hepdf->marginal(*mrv);
	shared_ptr<mpdf> c = hepdf->condition(crv);
	mepdf mm(m);

	Array<mpdf *> aa(2);
	aa(0) = c.get();
	aa(1) = &mm;
	mprod mEp(aa);

	mat smp = mEp.samplecond(vec(0), nsamples);
	vec emu = sum(smp, 2) / nsamples;

	// simplify overloading for Visual Studio
	vec delta = sqrt(variance) / sqrt(static_cast<double>(nsamples));
	CHECK_CLOSE_EX(mean, emu, delta);

	if (R.rows() > 0) {
	    mat er = (smp * smp.T()) / nsamples - outer_product(emu, emu);
	    CHECK_CLOSE_EX(R, er, tolerance);
	}

	// test of pdflog at zero
	vec zero(0);
	vec zeron(hepdf->dimension());
	for (int i = 0; i < zeron.size(); ++i) {
	    zeron(i) = 0;
	}

	CHECK_CLOSE_EX(hepdf->evallog(zeron), mEp.evallogcond(zeron, zero), tolerance);
    }
}

}