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

namespace bdm {

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

	UIFile in ( config_file_name );
	Array<shared_ptr<epdf_harness> > input;
	UI::get ( input, in, "data", UI::compulsory );
	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::optional );
	UI::get (nbins, set, "nbins", UI::optional );
	UI::get (nsamples, set, "nsamples", UI::optional );
	UI::get (R, set, "R", UI::optional );

	mrv = UI::build<RV> (set, "marginal_rv", UI::optional );

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

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 ) {
		int old_size = nbins.size();
		if ( old_size < 2 ) {
			ivec new_nbins ( "100 100" );
			for ( int i = 0; i < old_size; ++i ) {
				new_nbins ( i ) = nbins ( i );
			}

			nbins = new_nbins;
		}

		check_support_mean();
		check_support_integral();
	}

	if ( R.rows() > 0 ) {
		check_sample_mean();
		check_covariance();
	}

	if ( mrv ) {
		RV crv = hepdf->_rv().subt ( *mrv );
		epdf_ptr m = hepdf->marginal ( *mrv );
		mpdf_ptr c = hepdf->condition ( crv );

		mpdf_array aa ( 2 );
		aa ( 0 ) = c;
		aa ( 1 ) = new mepdf ( m );
		mprod mEp ( aa );

		check_cond_mean(mEp);
		if ( R.rows() > 0 ) {
			check_cond_covariance( mEp );
		}

		// test of pdflog at zero
		vec zero ( 0 );
		vec zeron=zeros ( hepdf->dimension() );

		double lpz = hepdf->evallog ( zeron );
		double lpzc=mEp.evallogcond ( zeron, zero );
		CHECK_CLOSE_EX ( lpz, lpzc, tolerance );

		vec lpzv(1);
		lpzv(0) = lpz;

		mat zero1n ( hepdf->dimension(), 1 );
		for ( int i = 0; i < zero1n.rows(); ++i ) {
			zero1n ( i, 0 ) = 0;
		}

		vec lpzv_act = hepdf->evallog_m ( zero1n );
		CHECK_CLOSE_EX ( lpzv, lpzv_act, tolerance );

		Array<vec> zeroa(3);
		lpzv = vec( zeroa.size() );
		for ( int i = 0; i < zeroa.size(); ++i ) {
			zeroa(i) = zeron;
			lpzv(i) = lpz;
		}

		lpzv_act = hepdf->evallog_m ( zeroa );
		CHECK_CLOSE_EX ( lpzv, lpzv_act, tolerance );
	}
}

void epdf_harness::check_support_mean() {
	vec xb = support.get_row ( 0 );
	vec yb = support.get_row ( 1 );

	int tc = 0;
	Array<vec> actual(CurrentContext::max_trial_count);
	do {
		vec emu = num_mean2 ( hepdf.get(), xb, yb, nbins ( 0 ), nbins ( 1 ) );
		actual( tc ) = emu;
		++tc;
	} while ( ( tc < CurrentContext::max_trial_count ) &&
		  !UnitTest::AreClose ( mean, actual( tc - 1 ), tolerance ) );
	if ( ( tc == CurrentContext::max_trial_count ) &&
	     ( !UnitTest::AreClose ( mean, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
		UnitTest::MemoryOutStream stream;
		stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << tolerance << " but was " << actual;

		UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);

		UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
       }
}

void epdf_harness::check_support_integral() {
	vec xb = support.get_row ( 0 );
	vec yb = support.get_row ( 1 );

	int tc = 0;
	Array<double> actual(CurrentContext::max_trial_count);
	do {
		double nc = normcoef ( hepdf.get(), xb, yb, nbins ( 0 ), nbins ( 1 ) );
		actual( tc ) = nc;
		++tc;
	} while ( ( tc < CurrentContext::max_trial_count ) &&
		  !UnitTest::AreClose ( 1.0, actual( tc - 1 ), tolerance ) );
	if ( ( tc == CurrentContext::max_trial_count ) &&
	     ( !UnitTest::AreClose ( 1.0, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
		UnitTest::MemoryOutStream stream;
		stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << tolerance << " but was " << actual;

		UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);

		UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
       }
}

void epdf_harness::check_sample_mean() {
	vec delta = make_close_tolerance ( variance, nsamples );

	int tc = 0;
	Array<vec> actual(CurrentContext::max_trial_count);
	do {
		mat smp = hepdf->sample_m ( nsamples );
		vec emu = smp * ones ( nsamples ) / nsamples;
		actual( tc ) = emu;
		++tc;
	} while ( ( tc < CurrentContext::max_trial_count ) &&
		  !UnitTest::AreClose ( mean, actual( tc - 1 ), delta ) );
	if ( ( tc == CurrentContext::max_trial_count ) &&
	     ( !UnitTest::AreClose ( mean, actual( CurrentContext::max_trial_count - 1 ), delta ) ) ) {
		UnitTest::MemoryOutStream stream;
		stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << delta << " but was " << actual;

		UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);

		UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
       }
}

void epdf_harness::check_covariance() {
	int tc = 0;
	Array<mat> actual(CurrentContext::max_trial_count);
	do {
		mat smp = hepdf->sample_m ( nsamples );
		vec emu = smp * ones ( nsamples ) / nsamples;
		mat er = ( smp * smp.T() ) / nsamples - outer_product ( emu, emu );
		actual( tc ) = er;
		++tc;
	} while ( ( tc < CurrentContext::max_trial_count ) &&
		  !UnitTest::AreClose ( R, actual( tc - 1 ), tolerance ) );
	
	if ( ( tc == CurrentContext::max_trial_count ) &&
	     ( !UnitTest::AreClose ( R, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
		UnitTest::MemoryOutStream stream;
		stream << CurrentContext::format_context(__LINE__) << "expected " << R << " +/- " << tolerance << " but was " << actual;

		UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);

		UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
       }
}

void epdf_harness::check_cond_mean( mprod &mep ) {
	vec delta = make_close_tolerance ( variance, nsamples );

	int tc = 0;
	Array<vec> actual(CurrentContext::max_trial_count);
	do {
		mat smp = mep.samplecond_m ( vec ( 0 ), nsamples );
		vec emu = sum ( smp, 2 ) / nsamples;
		actual( tc ) = emu;
		++tc;
	} while ( ( tc < CurrentContext::max_trial_count ) &&
		  !UnitTest::AreClose ( mean, actual( tc - 1 ), delta ) );
	if ( ( tc == CurrentContext::max_trial_count ) &&
	     ( !UnitTest::AreClose ( mean, actual( CurrentContext::max_trial_count - 1 ), delta ) ) ) {
		UnitTest::MemoryOutStream stream;
		stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << delta << " but was " << actual;

		UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);

		UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
       }
}

void epdf_harness::check_cond_covariance( mprod &mep ) {
	int tc = 0;
	Array<mat> actual(CurrentContext::max_trial_count);
	do {
		mat smp = mep.samplecond_m ( vec ( 0 ), nsamples );
		vec emu = sum ( smp, 2 ) / nsamples;
		mat er = ( smp * smp.T() ) / nsamples - outer_product ( emu, emu );
		actual( tc ) = er;
		++tc;
	} while ( ( tc < CurrentContext::max_trial_count ) &&
		  !UnitTest::AreClose ( R, actual( tc - 1 ), tolerance ) );
	if ( ( tc == CurrentContext::max_trial_count ) &&
	     ( !UnitTest::AreClose ( R, actual( CurrentContext::max_trial_count - 1 ), tolerance ) ) ) {
		UnitTest::MemoryOutStream stream;
		stream << CurrentContext::format_context(__LINE__) << "expected " << mean << " +/- " << tolerance << " but was " << actual;

		UnitTest::TestDetails details(*UnitTest::CurrentTest::Details(), 0, false);

		UnitTest::CurrentTest::Results()->OnTestFailure ( details, stream.GetText() );
       }
}

}