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

int main()
{
	using namespace heinlib;
	using namespace distributions;
	typedef double Real;

	srand(5);
	unsigned power_trials = 1000;

	static_matrix<3,1,Real> mean1 = {0,0,0};
	static_matrix<3,3,Real> covar1 = {1,0,0,
		0,1,0,
		0,0,1 };


	multivariate_normal<3,Real> X(mean1, covar1);


	unsigned num_samples = 10;


	std::ofstream shift_file("mv_mean_shift_effect.txt");
	xstream shift_out(shift_file);
	shift_out << "#mean_shift H_test_power E_test_power\n";

	//mean shift effect....
	std::cout << "mean shifting effect..\n";
	for(Real effect = 10; effect > 0.5; effect *= 0.8 )
	{
		//want the length of shift to be effect siuze...
		static_matrix<3,1,Real> mean2 = { effect/sqrt(3), effect/sqrt(3), effect/sqrt(3) };
		multivariate_normal<3,Real> Y(mean2, covar1);

		shift_out << effect << " ";

		unsigned H_reject = 0, E_reject=0, NN_reject=0, KS_reject=0;
		for(unsigned trial=0; trial<power_trials; ++trial)
		{
			std::vector< static_matrix<3,1,Real> > samplesX, samplesY;
			for(unsigned i=0; i<num_samples; ++i){
				samplesX.push_back(X.sample());
				samplesY.push_back(Y.sample());
			}

			Real p;
			p = bootstrap_pvalue( intersection_union_functor< multivariate_normal<3,Real> >(), samplesX, samplesY );
			if( p < 0.05 ) ++H_reject;

			p = bootstrap_pvalue( e_distance_functor(), samplesX, samplesY );
			if( p < 0.05 ) ++E_reject;

			p = bootstrap_pvalue( nearest_neighbor_functor(), samplesX, samplesY );
			if( p < 0.05 ) ++NN_reject;
		}

		shift_out << H_reject/(float)power_trials << " "
			      << E_reject/(float)power_trials << " "
				  << NN_reject/(float)power_trials << "\n";
	}


	std::ofstream scale_file("mv_covar_scale_effect.txt");
	xstream scale_out(scale_file);
	scale_out << "#covar_scale H_test_power E_test_power\n";

	//covar scale effect...
	std::cout << "covariance scaling effect..\n";
	for(Real effect = 10; effect > 0.5; effect *= 0.8 )
	{
		//want the length of shift to be effect siuze...
		static_matrix<3,3,Real> covar2 = covar1*effect;
		multivariate_normal<3,Real> Y(mean1, covar2);

		scale_out << effect << " ";

		unsigned H_reject = 0, E_reject=0, NN_reject=0, KS_reject=0;
		for(unsigned trial=0; trial<power_trials; ++trial)
		{
			std::vector< static_matrix<3,1,Real> > samplesX, samplesY;
			for(unsigned i=0; i<num_samples; ++i){
				samplesX.push_back(X.sample());
				samplesY.push_back(Y.sample());
			}

			Real p;
			p = bootstrap_pvalue( intersection_union_functor< multivariate_normal<3,Real> >(), samplesX, samplesY );
			if( p < 0.05 ) ++H_reject;

			p = bootstrap_pvalue( e_distance_functor(), samplesX, samplesY );
			if( p < 0.05 ) ++E_reject;

			p = bootstrap_pvalue( nearest_neighbor_functor(), samplesX, samplesY );
			if( p < 0.05 ) ++NN_reject;
		}

		scale_out << H_reject/(float)power_trials << " "
			      << E_reject/(float)power_trials << " "
				  << NN_reject/(float)power_trials << "\n";
	}


	/*
	std::ofstream rotate_file("mv_covar_rotate_effect.txt");
	xstream rotate_out(rotate_file);
	rotate_out << "#covar_rotate H_test_power E_test_power NN_test_power\n";

	//covar rotate effect...
	std::cout << "covariance rotate effect..\n";
	static_matrix<3,3,Real> covarR1 = {100,0,0,
		0,1,0,
		0,0,1 };
	multivariate_normal<3,Real> XR(mean1, covarR1);

	for(Real effect = 45; effect > 0.5; effect *= 0.8 )
	{
		//want the length of shift to be effect siuze...
		static_matrix<3,3,Real> R;
		static_matrix<3,1,Real> axis = {0,0,1}; //normalize( random_vector<3,Real>() );
		matrix_from_axis_angle( R, axis, effect*constants<Real>::deg_to_rad() );
		static_matrix<3,3,Real> covar2 = rotate_covariance( covarR1, R );

		//std::cout << "R = " << R << "\n";
		//std::cout << covar2 << "\n";

		multivariate_normal<3,Real> Y(mean1, covar2);

		rotate_out << effect << " ";

		unsigned H_reject = 0, E_reject=0, NN_reject=0, KS_reject=0;
		for(unsigned trial=0; trial<power_trials; ++trial)
		{
			std::vector< static_matrix<3,1,Real> > samplesX, samplesY;
			for(unsigned i=0; i<num_samples; ++i){
				samplesX.push_back(XR.sample());
				samplesY.push_back(Y.sample());
			}

			Real p;
			p = bootstrap_pvalue( intersection_union_functor< multivariate_normal<3,Real> >(), samplesX, samplesY );
			if( p < 0.05 ) ++H_reject;

			p = bootstrap_pvalue( e_distance_functor(), samplesX, samplesY );
			if( p < 0.05 ) ++E_reject;

			p = bootstrap_pvalue( nearest_neighbor_functor(), samplesX, samplesY );
			if( p < 0.05 ) ++NN_reject;
		}

		rotate_out << H_reject/(float)power_trials << " "
			      << E_reject/(float)power_trials << " "
				  << NN_reject/(float)power_trials << "\n";
	}*/
}