#include "test_util.h"
#include "stat/exp_family.h"
#include <fstream>

namespace bdm {

using itpp::vec;
using itpp::mat;

std::string load_test_file ( const char *fname ) {
	char buffer[8192];
	memset ( buffer, 0, sizeof ( buffer ) );
	std::ifstream src ( fname );
	src.read ( buffer, sizeof ( buffer ) - 1 );
	return std::string ( buffer );
}

double normcoef ( const epdf *ep, const vec &xb, const vec &yb,
                  int xn, int yn ) {
	mat Pdf ( xn + 1, yn + 1 );
	vec rgr ( 2 );

	double xstep = ( xb ( 1 ) - xb ( 0 ) ) / xn;
	double ystep = ( yb ( 1 ) - yb ( 0 ) ) / yn;

	double x = xb ( 0 );
	for ( int i = 0; i <= xn; x += xstep, i++ ) {
		rgr ( 0 ) = x;
		double y = yb ( 0 );
		for ( int j = 0; j <= yn; y += ystep, j++ ) {
			rgr ( 1 ) = y;
			Pdf ( i, j ) = exp ( ep->evallog ( rgr ) );
		}
	}

	return sumsum ( Pdf ) * xstep * ystep;
}

vec num_mean2 ( const epdf *ep, const vec &xb, const vec &yb,
                int xn, int yn ) {
	mat Pdf ( xn + 1, yn + 1 );
	vec rgr ( 2 );

	double xstep = ( xb ( 1 ) - xb ( 0 ) ) / xn;
	double ystep = ( yb ( 1 ) - yb ( 0 ) ) / yn;

	vec Mu ( xn + 1 );
	vec Si ( yn + 1 );

	double x = xb ( 0 );
	for ( int i = 0; i <= xn; x += xstep, i++ ) {
		Mu ( i ) = x;
		rgr ( 0 ) = x;
		double y = yb ( 0 );
		for ( int j = 0; j <= yn; y += ystep, j++ ) {
			Si ( j ) = y;
			rgr ( 1 ) = y;
			Pdf ( i, j ) = exp ( ep->evallog ( rgr ) );
		}
	}

	vec fm = sum ( Pdf, 2 );
	double sfm = sum ( fm );
	vec fs = sum ( Pdf, 1 );
	double sfs = sum ( fs );
	double vi0 = Mu * fm / sfm;
	double vi1 = Si * fs / sfs;
	return vec_2 ( vi0, vi1 );
}

}
