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

}