#include "gen_dirch.h"
#include <assert.h>
#include <math.h>
#include <time.h>
#include <unistd.h>	/* used for getpid to initialize random */
#include <sys/resource.h>	// for getrusage
#include <iostream>
using namespace std;

#include "Utilities/Random.h"
#include "gen_dirch.h"

double user_seconds(void)
{   struct rusage resources;
    getrusage(RUSAGE_SELF, &resources);
    return resources.ru_utime.tv_sec +
    	 resources.ru_utime.tv_usec*1.e-6;

}

inline static double abs(double x)
{    return x>=0? x: 0-x;
}

bool ok_dirch(const double* a, int numd, int c)
{
    double probs[numd];
    double sum[numd];
    double sum2[numd];
    
    double sum_a = 0.;
    for (int d=0; d<numd; d++)  
    {   sum_a += a[d]; 
        sum[d] = sum2[d] = 0.0;
    }
    gen_dirch g(numd, a);

    
    double start_time = user_seconds();
    for (int i=0; i<c; i++)
    {   g.generate(probs);
        for  (int d=numd-1; d>=0; d--)
	{   sum[d] += probs[d];
	    sum2[d] += probs[d]*probs[d];
	}       
    }
    double stop_time = user_seconds();

    bool ok=1;
    
    cout << "For sample of " << c << ", " << (stop_time-start_time)*1.e06/c
    	<< " microseconds/sample\n";
    for (int d=0; d<numd; d++)
    {	double mean= sum[d]/c;
        double e_mean=a[d]/sum_a;
	double err1 = mean-e_mean;
    	
        double square = sum2[d]/c;
	double e_square = a[d]*(a[d]+1)/ (sum_a*(sum_a+1));
	double err2 = square-e_square;
	
	cout << d << ": alpha=" << a[d] << " mean=" << mean
    		<< " expected_mean=" << e_mean 
		<< " err=" << err1
		<< "\n";
	cout << d << ":     mean square=" << square
		<< "  expected_mean_ square=" 
		<< e_square
		<< " err=" << err2
		<< "\n";
	if (abs(err1)>5.e-3 || abs(err2)>5.e-3)
	{   ok=0;
	    cout << "ERROR: error outside tolerance\n";
	}
    }
    cout << "\n" << flush;
    return ok;
}


int main()
{
    double beta[2] = {10., 100.};
    double tiny[20] = {0.153384, 0.0520756, 0.0073824, 0.0158439,
		0.428964, 0.025533, 0.0185789, 0.845361, 0.0282996,
		2.42256, 0.424296,  0.0190716, 0.0313429, 0.0274578,
		0.0252186, 0.028514, 0.0519217, 0.522946, 0.0279653,
		0.0664755};  
    double huge[20] = {64.0031, 15.1181, 59.4313, 67.6742, 47.888,
		66.0895, 26.101, 60.5549, 67.6252, 87.1146, 25.8355,
		46.3267, 47.446, 34.4312, 50.0487, 60.9615, 53.0045,
		66.1196, 16.4393, 41.099};  
    
    set_random(getpid());
    bool all_ok=1;
    
    all_ok &= ok_dirch(beta, 2, 2000000);
    all_ok &= ok_dirch(tiny, 20, 200000);
    all_ok &= ok_dirch(huge, 20, 200000);
    
    return all_ok? 0: 1;
}


// CHANGE LOG:
// Wed Jul 20 17:12:42 PDT 2005 Kevin Karplus
//	Chaged to C++
// Wed Jul 20 17:21:45 PDT 2005 Kevin Karplus
//	Added return value for all tests within tolerance
// Fri Dec 14 11:48:04 PST 2007 Kevin Karplus
//	added definition of abs