#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_mix.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;

}


/* This structure stores parameters of the Dirichlet mixture.
 * It is a stripped-down version of what is used in the SAM software.
 */
typedef struct{
  	int	AlphaChar;	/* Number of alphabet characters */
	int	L;		/* Number of prior distributions */
	double	*Mix;		/* Mixture coefficents for each prior */
	double	**Distr;	/* Prior distributions. L Dirchlet's 
				   over AlphaChar positions:
				   Distribution[L][AlphaChar+1] */
	} DirchMix;
	
void test_dirch_mix(const DirchMix *dm, int count)
{
    gen_dirch_mix g(dm->AlphaChar, dm->L, dm->Mix, 
    	const_cast<const double **>(dm->Distr));
    double probs[dm->AlphaChar];
    double sum[dm->AlphaChar];
    double sum2[dm->AlphaChar];
    for (int d=0; d<dm->AlphaChar; d++)
    {    sum[d] = sum2[d] = 0;
    }

    
    double sum_a[dm->L];
    
    double total_mix=0;
    for (int c=0; c< dm->L; c++)
    {   total_mix += dm->Mix[c];
	sum_a[c] = 0;
	for (int d=0; d<dm->AlphaChar; d++)
	{    sum_a[c] += dm->Distr[c][d+1];
	}
    }

    double start_time = user_seconds();;
    
    for (int i=0; i<count; i++)
    {   g.generate(probs);
        for  (int d=dm->AlphaChar-1; d>=0; d--)
	{   sum[d] += probs[d];
	    sum2[d] += probs[d]*probs[d];
	}       
    }
    double stop_time = user_seconds();
    
    cout << "For sample of " << count << ", "
    	<< (stop_time-start_time)*1.e06/count
	<< " microseconds/sample\n";
    for (int d=0; d<dm->AlphaChar; d++)
    {	double mixed_sum=0;
    	double mixed_sum2 = 0;
        for (int c=0; c< dm->L; c++)
	{   double a =  dm->Distr[c][d];
	    mixed_sum += dm->Mix[c]*a/sum_a[c];
	    mixed_sum2 += dm->Mix[c] * a*(a+1) / (sum_a[c]*(1+sum_a[c]));
	}

        cout << d << ": mean=" << sum[d]/count
		<< ", expected mean=" <<  mixed_sum/total_mix
		<< "\n";
    	cout << d << ": mean square=" << sum2[d]/count
		 << ", expected mean square="
		 <<  mixed_sum2/total_mix
		 << "\n";
    }
    cout << "\n";
}


int main()
{
    static double recode3_00[20]={0.668781, 0.072626, 1.08101, 0.930936,
	0.1235, 0.982387, 0.366701, 0.0586812, 0.849019, 0.185073,
	0.0648062, 0.979541, 0.519704, 0.555784, 0.636552, 1.04987,
	0.60931, 0.129148, 0.0676805, 0.201772}; 
    static double recode3_01[20]={0.880009, 0.197393, 0.303744, 0.497357,
	0.657492, 0.391988, 0.348159, 0.936094, 0.527206, 1.42844, 0.449302,
	0.368825, 0.384576, 0.439775, 0.581516, 0.622624, 0.747064, 1.08007,
	0.235012, 0.606928}; 
    static double recode3_02[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}; 
    static double recode3_03[20]={0.794007, 0.0130659, 0.624236, 1.85769,
	0.0290214, 0.115707, 0.123504, 0.22099, 1.52605, 0.341371, 0.111114,
	0.308302, 0.263545, 0.953727, 0.933444, 0.554741, 0.604551, 0.396451,
	0.00823516, 0.0420054}; 
    static double recode3_04[20]={0.740015, 0.187165, 0.0213261, 0.0456854,
	0.118944, 0.0633687, 0.0170331, 1.06684, 0.0380614, 0.733524,
	0.138456, 0.0300644, 0.0718692, 0.0240143, 0.0301022, 0.0862989,
	0.367283, 1.70735, 0.0113856, 0.045079}; 
    static double recode3_05[20]={0.15978, 0.0261585, 0.0505181, 0.125524,
	0.0350331, 0.102549, 0.157461, 0.0795041, 1.26261, 0.189383,
	0.0550608, 0.171028, 0.0844169, 0.290476, 1.44604, 0.129158, 0.138972,
	0.0851144, 0.0159134, 0.0637679}; 
    static double recode3_06[20]={0.308434, 0.0137217, 1.69731, 1.92422,
	0.0361113, 0.162357, 0.07232, 0.0487895, 0.236135, 0.0809074,
	0.0286236, 0.213663, 0.181631, 0.320245, 0.104878, 0.218398, 0.141668,
	0.0747719, 0.0141705, 0.0453433}; 
    static double recode3_07[20]={0.00260287, 9.99856e-06, 0.00631292,
	0.00445502, 0.00274753, 1.03886e-05, 1.02839e-05, 0.000913052,
	0.0029241, 0.00353485, 0.00105128, 0.00338172, 1.04172e-05,
	0.00173574, 0.00459583, 0.00274255, 0.00247625, 0.00175366,
	1.02411e-05, 0.00288489}; 
    static double recode3_08[20]={1.61043, 0.15522, 0.0378292, 0.0498243,
	0.0406484, 0.529136, 0.0217524, 0.040597, 0.0413396, 0.100193,
	0.0509779, 0.0357917, 0.0931204, 0.0367156, 0.0330646, 0.529587,
	0.196607, 0.230878, 0.00909518, 0.0329275}; 
    static double recode3_09[20]={0.15525, 0.0136827, 0.0857138, 0.0508316,
	0.0151451, 3.10555, 0.027169, 0.0140491, 0.0654038, 0.0257501,
	0.00901049, 0.127437, 0.0423873, 0.0345064, 0.0477247, 0.12452,
	0.0341196, 0.0230637, 0.00930115, 0.0187464}; 
    static double recode3_10[20]={0.225739, 0.0684326, 0.101072, 0.0813791,
	0.0298832, 0.0915218, 0.0336807, 0.0833114, 0.0931673, 0.0731542,
	0.0419314, 0.230216, 0.087446, 0.0694702, 0.0751969, 1.13857, 1.63158,
	0.179083, 0.00912576, 0.0311963}; 
    static double recode3_11[20]={1.3431e-06, 0.0166656, 0.00743068,
	1.34592e-06, 0.169076, 0.00407061, 0.00714122, 1.98221, 0.017522,
	0.816669, 0.114773, 0.00678027, 0.0106392, 0.0100244, 0.0158968,
	0.00879658, 0.0399043, 1.81043, 0.0150671, 0.0517801}; 
    static double recode3_12[20]={0.063525, 0.0288391, 1.09265, 0.0959581,
	0.00965196, 0.216914, 0.0730986, 0.0207325, 0.08719, 0.0315107,
	0.00960293, 0.752755, 0.059914, 0.0445321, 0.0312317, 0.287327,
	0.116896, 0.0249446, 0.00701663, 0.0305859}; 
    static double recode3_13[20]={0.294281, 0.019271, 0.12293, 0.162747,
	0.0373667, 0.145029, 0.0412349, 0.0815261, 0.157594, 0.151631,
	0.021412, 0.0601581, 3.6966, 0.0809085, 0.101856, 0.23533, 0.135424,
	0.140532, 0.00900473, 0.0321389}; 
    static double recode3_14[20]={0.0844832, 0.0584945, 0.0411628, 0.045719,
	0.847822, 0.0590839, 0.250253, 0.0675757, 0.0562614, 0.168617,
	0.0439737, 0.0794234, 0.028301, 0.0305672, 0.0598024, 0.0798202,
	0.0585385, 0.0858243, 0.227395, 1.30336}; 
    static double recode3_15[20]={0.0634034, 0.0246167, 1.3443e-06,
	0.00389272, 1.12953, 0.00796028, 1.35032e-06, 0.233395, 1.34466e-06,
	0.541933, 0.101309, 1.36412e-06, 0.027467, 0.00704479, 0.00802297,
	0.0248977, 0.0276933, 0.185467, 0.183309, 0.516892}; 
    static double recode3_16[20]={0.123696, 0.0454619, 0.0386434, 0.351847,
	0.0560181, 0.0439442, 0.223229, 0.01302, 0.148699, 0.19001, 0.120964,
	0.098734, 5.90055e-06, 0.554971, 0.219233, 0.0453885, 0.0564686,
	0.0614792, 0.0410248, 0.0800036}; 
    static double recode3_17[20]={0.0212037, 3.18769, 0.00745627, 0.00382411,
	0.00691924, 0.0126233, 1.34375e-06, 0.00724293, 0.00522979,
	0.00785563, 0.00489521, 0.0105326, 0.0136265, 0.00505819, 0.00677712,
	0.0251744, 0.0235516, 0.0371462, 0.00187667, 0.0038893}; 
    static double recode3_18[20]={0.0229376, 0.00427768, 0.00959934,
	0.013608, 0.182277, 0.0227654, 0.0157344, 0.0226783, 0.011561,
	0.0803491, 0.0154283, 0.00899225, 0.00980608, 0.00600945, 0.0342359,
	0.0216842, 0.0189306, 0.0223176, 1.83914, 0.154565}; 
    static double recode3_19[20]={2.16602e-06, 2.16245e-06, 0.0198496,
	2.17942e-06, 0.0246741, 2.47051e-06, 1.02563, 0.0131152, 2.16539e-06,
	0.00637704, 2.1414e-06, 0.0839371, 0.0168135, 0.0438887, 0.0252951,
	0.0235533, 0.0130626, 0.00797507, 2.16433e-06, 0.0545531}; 

    static double recode3_mix[20]= {0.176513, 0.207622, 0.0669246, 0.0868259,
    	0.0593123, 0.0358616, 0.03427, 0.0428319, 0.047875, 0.0466614,
    	0.0283695, 0.0301127, 0.0233828, 0.034662, 0.0270202, 0.0226822,
    	0.00898452, 0.00716226, 0.00710292, 0.00582299}; 
    
    static double *recode3_comps[20] = {
    	recode3_00, recode3_01, recode3_02, recode3_03, recode3_04,
    	recode3_05, recode3_06, recode3_07, recode3_08, recode3_09, 
    	recode3_10, recode3_11, recode3_12, recode3_13, recode3_14,
    	recode3_15, recode3_16, recode3_17, recode3_18, recode3_19};
    DirchMix recode3 = {20, 20, recode3_mix, recode3_comps};
    	    
    
    
    static double mix_coeff[6] = {0.272155, 0.259417, 0.170403, 0.139492,
			    0.0813267, 0.0772063}; 

    static double comp0[20]={ 44.1267, 8.54778, 34.3287, 45.012, 25.4885,
	    40.1706, 13.8165, 37.9695, 40.2287, 57.0279, 15.7816, 24.6325,
	    26.6716, 21.0461, 33.2734, 39.015, 31.1269, 43.4722, 6.95188,
	    19.8187}; 
    static double comp1[20]={ 12.9502, 3.56352, 9.54209, 10.7789,
	    7.22886, 11.7739, 4.45173, 9.25549, 10.113, 16.1048,
	    4.26184, 8.89566, 
	    10.1378, 8.39929, 9.83732, 16.2176, 11.391, 11.4088, 2.43706,
	    5.69262}; 
    static double comp2[20]={ 17.2777, 5.02666, 20.633, 27.1002, 18.6022,
	    16.5669, 7.46735, 29.459, 33.201, 36.5843, 8.33165,
	    23.2366, 12.2568, 
	    12.8126, 15.8195, 26.7285, 18.5344, 20.7154, 3.5349, 14.591}; 
    static double comp3[20]={ 36.82, 3.69397, 16.532, 17.1343, 9.09368,
	    25.9803, 7.44929, 11.6098, 7.16141, 31.5673, 6.2929,
	    6.50351, 17.932, 
	    9.98346, 24.4647, 17.2535, 16.1794, 23.996, 4.379, 6.65627}; 
    static double comp4[20]={ 2.4706, 0.588751, 1.33572, 1.65825,
	    1.08549, 2.23149, 0.722458, 1.4641, 1.70528, 2.49194, 0.900071,
	    1.33377, 1.84355, 1.25755, 1.65407, 2.63801, 1.88732, 1.84009,
	    0.399507, 0.92818}; 
    static double comp5[20]={ 25.0801, 4.5497, 8.56496, 10.2153, 24.5674,
	    21.9557, 5.42381, 29.6721, 13.1884, 43.4575, 10.6261, 11.4957,
	    13.0346, 8.10884, 11.8077, 25.3051, 18.1739, 25.4028, 5.88753,
	    12.6309}; 

    static double *comps[6] = {comp0, comp1, comp2, comp3, comp4, comp5};
    DirchMix rsdb_comp = {20,6, mix_coeff, comps};
    
    set_random(getpid());
    test_dirch_mix(&recode3, 200000);
    test_dirch_mix(&rsdb_comp, 200000);
    
    return 0;
}

// Thu Jul 21 08:48:25 PDT 2005 Kevin Karplus
//	No longer use awkward indexing starting at 1 for components.