#include "heinlibs.h"
#include <fstream>

int main()
{
	using namespace heinlib;
	using namespace distributions;

	typedef double Real;

	unsigned num_samples = 100;

	srand(5);

	normal<Real> X(1,0);
	std::vector<Real> samplesX;
	for(unsigned i=0; i<num_samples; ++i)
		samplesX.push_back(X.sample());

	/*
	for(Real m=-2; m<2; m += 0.1 )
	{
		normal<Real> Y(1,m);
		std::vector<Real> samplesY;
		for(unsigned i=0; i<num_samples; ++i)
			samplesY.push_back(Y.sample());

		std::cout << X << " vs. " << Y << " : H-test = " << bootstrap_pvalue( normal_equality( samplesX, samplesY ), samplesX, samplesY, false ) << "\n";
		//std::cout << X << " vs. " << Y << " : Q-Dist = " << bootstrap_pvalue( QDA_statistic( samplesX, samplesY ), samplesX, samplesY ) << "\n";
		//std::cout << X << " vs. " << Y << " : Q-Dist = " << QDA_statistic( samplesX, samplesY ) << "\n";
		//std::cout << X << " vs. " << Y << " : e-Dist = " << bootstrap_pvalue( e_distance( samplesX, samplesY ), samplesX, samplesY ) << "\n\n";
		//std::cout << X << " vs. " << Y << " : NN-dist = " << bootstrap_pvalue( nearest_neighbor_statistic( samplesX, samplesY ), samplesX, samplesY ) << "\n\n";
		//system("pause");
	}*/

	std::ofstream H_dist_file("H_contour.txt"), E_dist_file("E_contour.txt"), NN_dist_file("NN_contour.txt"), KS_contour("KS_contour.txt");
	xstream H_out(H_dist_file), E_out(E_dist_file), NN_out(NN_dist_file), KS_out(KS_contour);

	H_out << "# samples = " << num_samples << "\n# mean sigma stat\n";
	E_out << "# samples = " << num_samples << "\n# mean sigma stat\n";
	NN_out << "# samples = " << num_samples << "\n# mean sigma stat\n";
	KS_out << "# samples = " << num_samples << "\n# mean sigma stat\n";


	std::ofstream H_dist_filep("H_contourp.txt"), E_dist_filep("E_contourp.txt"), NN_dist_filep("NN_contourp.txt"), KS_contourp("KS_contourp.txt");
	xstream H_outp(H_dist_filep), E_outp(E_dist_filep), NN_outp(NN_dist_filep), KS_outp(KS_contourp);

	H_outp << "# samples = " << num_samples << "\n# mean sigma p-value\n";
	E_outp << "# samples = " << num_samples << "\n# mean sigma p-value\n";
	NN_outp << "# samples = " << num_samples << "\n# mean sigma p-value\n";
	KS_outp << "# samples = " << num_samples << "\n# mean sigma p-value\n";

	std::ofstream NPH_dist_file("NPH_contour.txt");
	xstream NPH_out(NPH_dist_file);

	NPH_out << "# samples = " << num_samples << "\n# mean sigma stat\n";

	std::ofstream NPH_dist_filep("NPH_contourp.txt");
	xstream NPH_outp(NPH_dist_filep);

	NPH_outp << "# samples = " << num_samples << "\n# mean sigma p-value\n";


	for(unsigned j=0; j<=15; ++j)
	{
		//0.1 to 1 in 5 steps, but we go beyond that to get up to approx s=10..
		Real a = 0.1 - 1;
		Real b = pow( 1+1-0.1, 1/5.0 );
		Real s = a + pow(b,(double)j);


		for(Real m=-2; m<2; m += 0.1 )
		{
			normal<Real> Y(s,m);
			std::vector<Real> samplesY;
			for(unsigned i=0; i<num_samples; ++i)
				samplesY.push_back(Y.sample());


			H_out << m << " " << s << " " << intersection_union_statistic< normal<Real> >( samplesX, samplesY ) << "\n";
			E_out << m << " " << s << " " << e_distance( samplesX, samplesY ) << "\n";
			NN_out << m << " " << s << " " << nearest_neighbor_statistic( samplesX, samplesY ) << "\n";
			KS_out << m << " " << s << " " << KS_statistic( samplesX, samplesY ) << "\n";
			NPH_out << m << " " << s << " " << intersection_union_statistic< kernel_density<1,Real> >( samplesX, samplesY ) << "\n";

			NPH_outp << m << " " << s << " " << bootstrap_pvalue( intersection_union_functor< kernel_density<1,Real> >(), samplesX, samplesY ) << "\n";
			H_outp << m << " " << s << " " << bootstrap_pvalue( intersection_union_functor< normal<Real> >(), samplesX, samplesY ) << "\n";
			E_outp << m << " " << s << " " << bootstrap_pvalue( e_distance_functor(), samplesX, samplesY ) << "\n";
			NN_outp << m << " " << s << " " << bootstrap_pvalue( nearest_neighbor_functor(), samplesX, samplesY ) << "\n";
			KS_outp << m << " " << s << " " << bootstrap_pvalue( twosample_KS_functor<Real>(), samplesX, samplesY ) << "\n";



		}
		H_out << "\n";
		E_out << "\n";
		NN_out << "\n";
		KS_out << "\n";

		H_outp << "\n";
		E_outp << "\n";
		NN_outp << "\n";
		KS_outp << "\n";

		NPH_out << "\n";
		NPH_outp << "\n";


	}
}