#include "gen_dirch.h"
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <time.h>
#include <unistd.h>	/* used for getpid to initialize random */

void test_dirch(const double* a, int numd, int c)
{
    int i, d;
    clock_t start_time,     stop_time;
    
    double *probs;
    double *sum, *sum2;
    double sum_a;
    
    gen_dirch_param g;
    
    gen_dirch_initialize(&g, numd, a);
    probs = (double*)calloc(numd, sizeof(double));
    sum = (double*) calloc(numd, sizeof(double));
    sum2= (double*) calloc(numd, sizeof(double));
    
    sum_a = 0.;
    for (d=0; d<numd; d++)
    {    sum_a += a[d];
    }

    start_time = clock();
    
    for (i=0; i<c; i++)
    {   gen_dirch(&g, probs);
        for  (d=numd-1; d>=0; d--)
	{   sum[d] += probs[d];
	    sum2[d] += probs[d]*probs[d];
	}       
    }
    stop_time = clock();
    
    printf("For sample of %g, %g microseconds/sample\n",	
    	(double)c, (stop_time-start_time+0.0)/c);
    for (d=0; d<numd; d++)
    {	printf("%d: alpha=%g mean= %g, expected mean=%g\n", 
    		d, a[d], sum[d]/c, a[d]/sum_a);
    	printf("%d: mean square= %g, expected mean square=%g\n", 
    		d, sum2[d]/c, a[d]*(a[d]+1)/ (sum_a*(sum_a+1)));
    }
    printf("\n");
}


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};  
    
    srandom(getpid());
    test_dirch(beta, 2, 2000000);
    test_dirch(tiny, 20, 200000);
    test_dirch(huge, 20, 200000);
    
}