/*
 * GMM.cl
 *
 *  Created on: Jan 7, 2012
 *      Author: scannon
 */
 
#ifndef GMM_CL
#define GMM_CL

//#include "ranluxcl.cl"
//#include "Perceptron.cl"
 
struct __GMM__;
typedef struct __GMM__
{
	// 1 input per conditional variable
	// 3 outputs per normal:  mean, standard deviation, and weight
	Perceptron *perceptron;
	struct __GMM__ * conditionals[16];
	int numConditionals;
	char sampled;
	float value;
} GMM;
 
 
GMM * GMM__init(GMM *self,
		Perceptron *perceptron,
		GMM * conditionals[16],
		int numConditionals)
{
	int i = 0;
	
	if(numConditionals > 16 || numConditionals < 0)
		return 0;
	
	self->perceptron = perceptron;
	for(i = 0; i < numConditionals; i++)
		self->conditionals[i] = conditionals[i];
	self->numConditionals = numConditionals;
	self->sampled = 0;
	
	return self;
}
 
/**
 * O(O(Perceptron_calc()) + 16 + 16 + O(ranluxcl32norm())) => O(1)
 * M(sizeof(Perceptron) + sizeof(ptr)*16 + sizeof(int) + sizeof(char) + sizeof(float)) => M(4k + 77)
 */
float GMM_sample(GMM *self, ranluxcl_state_t *rst)
{
	float normParams[48];
	int outNodeIndex = 0;
	float4 ranluxRet;
	float mean = 0.0f;
	float sd = 0.0f;
	float weightSum = 0.0f;
	float16 cVals;
	int i = 0;
	
	// Return if already sampled
	if(self->sampled)
		return self->value;
	
	// Sample parents first
	cVals.s0 = self->numConditionals >= 1 ? GMM_sample(self->conditionals[0], rst) : 0.0f;
	cVals.s1 = self->numConditionals >= 2 ? GMM_sample(self->conditionals[1], rst) : 0.0f;
	cVals.s2 = self->numConditionals >= 3 ? GMM_sample(self->conditionals[2], rst) : 0.0f;
	cVals.s3 = self->numConditionals >= 4 ? GMM_sample(self->conditionals[3], rst) : 0.0f;
	cVals.s4 = self->numConditionals >= 5 ? GMM_sample(self->conditionals[4], rst) : 0.0f;
	cVals.s5 = self->numConditionals >= 6 ? GMM_sample(self->conditionals[5], rst) : 0.0f;
	cVals.s6 = self->numConditionals >= 7 ? GMM_sample(self->conditionals[6], rst) : 0.0f;
	cVals.s7 = self->numConditionals >= 8 ? GMM_sample(self->conditionals[7], rst) : 0.0f;
	cVals.s8 = self->numConditionals >= 9 ? GMM_sample(self->conditionals[8], rst) : 0.0f;
	cVals.s9 = self->numConditionals >= 10 ? GMM_sample(self->conditionals[9], rst) : 0.0f;
	cVals.sA = self->numConditionals >= 11 ? GMM_sample(self->conditionals[10], rst) : 0.0f;
	cVals.sB = self->numConditionals >= 12 ? GMM_sample(self->conditionals[11], rst) : 0.0f;
	cVals.sC = self->numConditionals >= 13 ? GMM_sample(self->conditionals[12], rst) : 0.0f;
	cVals.sD = self->numConditionals >= 14 ? GMM_sample(self->conditionals[13], rst) : 0.0f;
	cVals.sE = self->numConditionals >= 15 ? GMM_sample(self->conditionals[14], rst) : 0.0f;
	cVals.sF = self->numConditionals >= 16 ? GMM_sample(self->conditionals[15], rst) : 0.0f;
 
	// Determine weights, means, and standard deviations from conditionals
	Perceptron_calc2(self->perceptron, cVals, normParams);
 
	// Determine sum of weight ratios
	weightSum = 0.0f;
	for(i = 0; i < 16; i++)
		weightSum += normParams[i*3 + 2];
 
	// Randomly choose normal from uniform distribution (0 to 1 exclusive)
	ranluxRet = ranluxcl32(rst);
	outNodeIndex = 0;
	while(ranluxRet.x > 0)
	{
		ranluxRet.x -= normParams[outNodeIndex*3 + 2]/weightSum;
		outNodeIndex++;
	}
	outNodeIndex--;
 
	// Sample normal distributions
	mean = normParams[3*outNodeIndex + 0];
	sd = normParams[3*outNodeIndex + 1];
	self->value = ranluxcl32norm(rst).x*sd + mean;
	self->sampled = 1;
	
	return self->value;
}

#endif	// GMM_CL