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

}