Monte Carlo Multilayer Perceptron Gaussian Mixture Model Distribution

The Monte Carlo Multilayer Perceptron Gaussian Mixture Model Distribution (MCMPGMMDist).

A Bayesian network node conditional probability distribution is represented by a series of Gaussian mixture models. A Gaussian mixture model is defined for every combination of states the parent nodes of the given node can take on. The parameters for each Gaussian mixture model is determined from a multilayer perceptron. As the mixture models for each set of parent states are not independent of each other, the neural network serves to also define and maintain the relationship of the mixture models of the conditional probability distribution.

Ultimately, the states of the parents of the given node are fed into the neural network. There is an input node for each parent node. The output nodes of the neural network then provide the parameters of the Gaussian mixture model for the given state set of the parent nodes. These parameters are the weight, mean, and standard deviation of each normal distribution in the Gaussian mixture model. This set of parameters define the Gaussian mixture model for the particular state set of the parent nodes of the given node.

The advantage of this methodology is that any continuous probability distribution can be approximated by the combination of the multilayer perceptron and the Gaussian mixture model. Furthermore, the relationship between the distributions for different parent state sets is maintained by the neural network. By this methodology recognizes this relationship, we are able to train this distribution more easily and represent and represent the nodes conditional probability distribution more compactly.

XOR Problem

x[1]

x[1, {'y':2, 'z':3}]

MCMPGMMDistTest.py

'''
Created on Nov 19, 2009
 
@author: Stephen
'''
import unittest
from MCMPGMMDist import MCMPGMMDist
from numpy import sqrt;
 
 
class Test(unittest.TestCase):
	def testAddConditional(self):
		input1 = [2, 3, 7];
		input2 = [2, 3, 7, 0];
		dist = MCMPGMMDist(['a', 'b'], 0, 10);
		output1 = dist.net(input1);
		dist.addConditional('c');
		output2 = dist.net(input2);
		self.assertEqual(output1, output2);
		self.assertEquals(dist.condVars, ['a', 'b', 'c']);
 
	def testRemoveConditional(self):
		input1 = [2, 3, 0, 7];
		input2 = [2, 3, 7];
		dist = MCMPGMMDist(['a', 'b', 'c'], 0, 10);
		output1 = dist.net(input1);
		dist.removeConditional('b');
		output2 = dist.net(input2);
		self.assertEqual(output1, output2);
		self.assertEquals(dist.condVars, ['a', 'c']);
 
	def testNumWeights(self):
		dist = MCMPGMMDist(['a', 'b', 'c'], 0, 10);
		self.assertEquals(dist.numWeights(), 5 * MCMPGMMDist.NUM_HIDDEN_VARS + (MCMPGMMDist.NUM_HIDDEN_VARS + 1) * MCMPGMMDist.NUM_NORMALS * 3);
 
	def testModifyWeightWithinBounds(self):
		# Setup test
		dist = MCMPGMMDist(['a', 'b', 'c'], 0, 10);
		wStartValue = dist.net.weights[1];
		bound = 2.38 / sqrt(5);
 
		# Modify weights a little and test
		dist.modifyWeightWithinBounds(1, 0.1);
		self.assertAlmostEquals(dist.net.weights[1] - wStartValue, (bound - wStartValue) * .1);
 
		# Test extremes of modifier
		dist.modifyWeightWithinBounds(1, 1);
		self.assertAlmostEquals(dist.net.weights[1], bound);
		dist.modifyWeightWithinBounds(1, -1);
		self.assertAlmostEquals(dist.net.weights[1], -bound);
		dist.modifyWeightWithinBounds(1, 0.5);
		self.assertAlmostEquals(dist.net.weights[1], 0);
 
		# Test exception thrown when going out of bounds
		self.assertRaises(ValueError, dist.modifyWeightWithinBounds, 1, 1.1);
 
	def testModifyWeight(self):
		dist = MCMPGMMDist(['a', 'b', 'c'], 0, 10);
		wStartValue = dist.net.weights[1];
		dist.modifyWeight(1, 4.4);
		self.assertAlmostEquals(dist.net.weights[1] - wStartValue, wStartValue * 4.4);
 
 
if __name__ == "__main__":
	#import sys;sys.argv = ['', 'Test.testAddConditional', 'Test.testRemoveConditional', 'Test.testNumWeights', 'Test.testModifyWeightWithinBounds', 'Test.testModifyWeight']
	unittest.main()

MCMPGMMDist.py

'''
Created on Aug 1, 2009
 
@author: Stephen
'''
 
from ffnet import ffnet, mlgraph, savenet, loadnet, exportnet;
from numpy.random import normal, random_sample;
from numpy import array, digitize, zeros, where, insert, delete, sqrt;
from copy import deepcopy
 
 
class MCMPGMMDist(object):
	'''
	Monte Carlo Multilayer Perceptron Gaussian Mixture Model Distribution
	'''
 
	NUM_HIDDEN_VARS = 10;
	NUM_NORMALS = 10;
 
	def __init__(self, condVars, normMin, normMax):
		'''
		Constructor
		condVars - List of identifiers for conditional variables for this conditional distribution
		normMin - Minimum possible value any of the conditional variables can be
		normMax - Maximum possible value any of the conditional variables can be
		'''
		if(condVars.__class__ != [].__class__):
			raise TypeError("concVars must be a list");
		self.condVars = condVars;
		self.numCondVars = len(self.condVars);
		self.normMin = normMin;
		self.normMax = normMax;
		conec = mlgraph((self.numCondVars + 1, MCMPGMMDist.NUM_HIDDEN_VARS, MCMPGMMDist.NUM_NORMALS * 3));
		self.net = ffnet(conec);
 
	def sample(self, *args, **keywords):
		'''
		Produces a sample from this distribution
		*args - Value of each of the variables this distribution is conditioned
		upon and x.  Must have a value for each conditional variable and x.
		numSamples - Number of samples to generate.  Default is 1000.
		'''
		# Ensure each conditional variable is provided
		assert(len(args) == self.numCondVars + 1);
 
		# Get number of samples if provided
		numSamples = 1000;
		if 'numSamples' in keywords:
			numSamples = keywords['numSamples'];
 
		# Normalize values of variables this distribution is conditioned upon
		inputs = [args[0]];
		if len(args) > 1:
			for condVals in args[1].keys():
				if condVals not in self.condVars:
					raise Exception("Conditional value not a parent of node.");
		for cond in self.condVars:
			inputs.append(args[1][cond]);
		values = (array(inputs) - self.normMin) / (self.normMax - self.normMin);
 
		# Get parameters for each normal distribution
		normParams = self.net(values);		# Need exact number of inputs in values
		normParams = [normParams[3 * i : 3 * i + 3] for i in range(len(normParams) / 3)];
		normParams = array(normParams).T;
		weights = normParams[0];
		weights = weights / sum(weights);
		normParams[0] = weights;
		normParams = normParams.T;
 
		# Generate sample
		retVal = [];
		cWeights = [sum(weights[0 : i + 1]) for i in range(len(weights))];
		bins = digitize(random_sample(numSamples), cWeights);
		mask = zeros([MCMPGMMDist.NUM_NORMALS, numSamples]);
		for i in range(numSamples):
			mask[bins[i]][i] = 1;
		for i in range(MCMPGMMDist.NUM_NORMALS):
			mean = normParams[i][1] * (self.normMax - self.normMin) + self.normMin;
			sd = normParams[i][2] * (self.normMax - self.normMin) + self.normMin;
			retVal += [normal(mean, sd, numSamples).T];
		retVal = retVal * mask;
 
		return retVal[where(retVal != 0)];
 
	def addConditional(self, newCondVar):
		'''
		Adds a conditional variable to this distribution
		newCondVar - Identifier for new conditional variable
		'''
		# net1 = ffnet(mlgraph((2, 3, 3), False))
		# net2 = ffnet(mlgraph((3, 3, 3), False))
		# net1.weights = array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])
		# net2.weights = array([1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])
		# net1([1, 1])
		# -> [0.97932062734407388, 0.97932062734407388, 0.97932062734407388]
		# net2([1, 1, 1])
		# -> [0.97932062734407388, 0.97932062734407388, 0.97932062734407388]
 
		# Create new neural network
		conec = mlgraph((self.numCondVars + 2, MCMPGMMDist.NUM_HIDDEN_VARS, MCMPGMMDist.NUM_NORMALS * 3));
		newNet = ffnet(conec);
 
		# Insert weights for new variable
		newNet.weights = self.net.weights;
		insertPoints = range(1, MCMPGMMDist.NUM_HIDDEN_VARS * (self.numCondVars + 2) + 1, self.numCondVars + 2);
		bound = 2.38 / sqrt(self.numCondVars + 3);
		newNet.weights = insert(newNet.weights, insertPoints, random_sample(len(insertPoints)) * bound * 2 - bound);
		self.net = newNet;
 
		# Update remaining variables
		self.condVars.append(newCondVar);
		self.numCondVars = len(self.condVars);
 
	def removeConditional(self, condVar):
		'''
		Removes a conditional variable from this distribution
		condVar - Identifier for the conditional variable to remove
		'''
		# Create new neural network
		conec = mlgraph((self.numCondVars, MCMPGMMDist.NUM_HIDDEN_VARS, MCMPGMMDist.NUM_NORMALS * 3));
		newNet = ffnet(conec);
 
		# Insert weights for new variable
		newNet.weights = self.net.weights;
		removalPoints = range((self.numCondVars - self.condVars.index(condVar) - 1) + 1, MCMPGMMDist.NUM_HIDDEN_VARS * (self.numCondVars + 2), self.numCondVars + 2);
		newNet.weights = delete(newNet.weights, removalPoints);
		self.net = newNet;
 
		# Update remaining variables
		self.condVars.remove(condVar);
		self.numCondVars = len(self.condVars);
 
	def numWeights(self):
		'''
		Returns the number of weights in the neural network backing this distribution
		'''
		return len(self.net.weights);
 
	def modifyWeightWithinBounds(self, wIndex, modifier):
		'''
		Modifies one of the weights of the neural network backing this distribution within
		certain bounds for those weights.  Due to the behavior of this method, the amount in which
		a unit modifier modifies the weight is a function of how far that weight is from its
		bounds.  Thus, a modifier of -0.1 will modify a positive weight more than a modifier of 0.1
		since a positive weight is farther from its lower bound.  In essence, a modifier of -0.1
		will move a weight 10% towards its lower bound while a modifier of 0.1 will move a weight
		10% towards its upper bound.
		wIndex - The index of the weight to modify.  Must be between 0 and the number of weights in
		the neural network backing this distribution.
		modifier - A value between -1 and 1 representing the amount to modify the weight towards
		its bounds.  A positive value will increase the weight while a negative value will decrease
		the weight.  A value of 1 will put the weight at its maximum bound while a value of -1 will
		put the weight at its minimum bound.  A value of 0 will not modify the weight and any value
		in between 0 and 1 or 0 and -1 will modify the value linearly towards its bound.  Values
		greater than 1 or -1 are not allowed.
		'''
		# Check bounding conditions of modifier
		if modifier > 1:
			raise ValueError, 'Modifier value greater than 1:  modifier = ' + str(modifier);
		if modifier < -1:
			raise ValueError, 'Modifier value less than -1:  modifier = ' + str(modifier);
 
		# Modify weight
		inputs = self.numCondVars + 2;
		bound = 0;
		if wIndex < MCMPGMMDist.NUM_HIDDEN_VARS * inputs:
			bound = 2.38 / sqrt(inputs);
		else:
			bound = 2.38 / sqrt(MCMPGMMDist.NUM_HIDDEN_VARS + 1);
		if modifier > 0:
			range = bound - self.net.weights[wIndex];
			self.net.weights[wIndex] += range * modifier;
		else:
			range = bound + self.net.weights[wIndex];
			self.net.weights[wIndex] -= range * -1 * modifier;
 
	def modifyWeight(self, wIndex, modifier):
		'''
		Modifies one of the weights of the neural network backing this distribution within
		certain bounds for those weights.
		wIndex - The index of the weight to modify.  Must be between 0 and the number of weights in
		the neural network backing this distribution.
		modifier - The modifier for the weight.  This value is the percent change to apply to the
		weight where 1 is 100% and -1 is -100%
		'''
		self.net.weights[wIndex] += self.net.weights[wIndex] * modifier;
 
	def __deepcopy__(self, memo):
		cpy = MCMPGMMDist(self.condVars[:], self.normMin, self.normMax);
		#cpy.condVars = self.condVars[:];
		#cpy.numCondVars = self.numCondVars;
		#cpy.normMin = self.normMin;
		#cpy.normMax = self.normMax;
		cpy.net = deepcopy(self.net, memo);
		return cpy;

Perceptron.cl

/*
 * Perceptron.cl
 *
 *  Created on: Jan 22, 2012
 *      Author: scannon
 */
 
 
typedef struct __PerceptronNode__
{
	float16 w;
} PerceptronNode;
 
typedef struct __Perceptron__
{
	int numInputNodes;
	PerceptronNode hiddenNodes[16];
	PerceptronNode outputNodes[48];
} Perceptron;
 
/**
 * O(1)
 * horizontal asymptotes at 0 and 1
 * sigmoid(0) = 0.5
 */
float sigmoid(float x)
{
	return 1.0f / (1.0f + exp(-1.0f * x));
}
 
PerceptronNode * PerceptronNode_init(PerceptronNode *self, float16 w)
{
	self->w = w;
 
	return self;
}
 
PerceptronNode * PerceptronNode_init2(PerceptronNode *self, int numInputNodes)
{
	self->w.s0 = numInputNodes <= 1 ? 1 : 0;
	self->w.s1 = numInputNodes <= 2 ? 1 : 0;
	self->w.s2 = numInputNodes <= 3 ? 1 : 0;
	self->w.s3 = numInputNodes <= 4 ? 1 : 0;
	self->w.s4 = numInputNodes <= 5 ? 1 : 0;
	self->w.s5 = numInputNodes <= 6 ? 1 : 0;
	self->w.s6 = numInputNodes <= 7 ? 1 : 0;
	self->w.s7 = numInputNodes <= 8 ? 1 : 0;
	self->w.s8 = numInputNodes <= 9 ? 1 : 0;
	self->w.s9 = numInputNodes <= 10 ? 1 : 0;
	self->w.sA = numInputNodes <= 11 ? 1 : 0;
	self->w.sB = numInputNodes <= 12 ? 1 : 0;
	self->w.sC = numInputNodes <= 13 ? 1 : 0;
	self->w.sD = numInputNodes <= 14 ? 1 : 0;
	self->w.sE = numInputNodes <= 15 ? 1 : 0;
	self->w.sF = numInputNodes <= 16 ? 1 : 0;
 
	return self;
}
 
PerceptronNode * PerceptronNode_init(PerceptronNode *self,
		int numInputNodes,
		ranluxcl_state_t *rst)
{
	self->w.s0 = numInputNodes <= 1 ? ranluxcl32(rst) * 20 - 10 : 0;
	self->w.s1 = numInputNodes <= 2 ? ranluxcl32(rst) * 20 - 10 : 0;
	self->w.s2 = numInputNodes <= 3 ? ranluxcl32(rst) * 20 - 10 : 0;
	self->w.s3 = numInputNodes <= 4 ? ranluxcl32(rst) * 20 - 10 : 0;
	self->w.s4 = numInputNodes <= 5 ? ranluxcl32(rst)  * 20 - 10: 0;
	self->w.s5 = numInputNodes <= 6 ? ranluxcl32(rst) * 20 - 10 : 0;
	self->w.s6 = numInputNodes <= 7 ? ranluxcl32(rst) * 20 - 10 : 0;
	self->w.s7 = numInputNodes <= 8 ? ranluxcl32(rst) * 20 - 10 : 0;
	self->w.s8 = numInputNodes <= 9 ? ranluxcl32(rst) * 20 - 10 : 0;
	self->w.s9 = numInputNodes <= 10 ? ranluxcl32(rst) * 20 - 10 : 0;
	self->w.sA = numInputNodes <= 11 ? ranluxcl32(rst) * 20 - 10 : 0;
	self->w.sB = numInputNodes <= 12 ? ranluxcl32(rst) * 20 - 10 : 0;
	self->w.sC = numInputNodes <= 13 ? ranluxcl32(rst) * 20 - 10 : 0;
	self->w.sD = numInputNodes <= 14 ? ranluxcl32(rst) * 20 - 10 : 0;
	self->w.sE = numInputNodes <= 15 ? ranluxcl32(rst) * 20 - 10 : 0;
	self->w.sF = numInputNodes <= 16 ? ranluxcl32(rst) * 20 - 10 : 0;
 
	return self;
}
 
/**
 * f(x) = sigmoid(sum(w[0]*x[0]...w[n-1]*x[n-1]))
 * f(x) = sigmoid(dot(w, x))
 * O(numInputs)
 */
float PerceptronNode_calc(PerceptronNode *self, float16 x)
{
	float sum = 0.0f;
 
	sum += self->w.s0*x.s0;
	sum += self->w.s1*x.s1;
	sum += self->w.s2*x.s2;
	sum += self->w.s3*x.s3;
	sum += self->w.s4*x.s4;
	sum += self->w.s5*x.s5;
	sum += self->w.s6*x.s6;
	sum += self->w.s7*x.s7;
	sum += self->w.s8*x.s8;
	sum += self->w.s9*x.s9;
	sum += self->w.sA*x.sA;
	sum += self->w.sB*x.sB;
	sum += self->w.sC*x.sC;
	sum += self->w.sD*x.sD;
	sum += self->w.sE*x.sE;
	sum += self->w.sF*x.sF;
 
	return sigmoid(sum);
}
 
void Perceptron_init(Perceptron *self,
		int numInputNodes,
		PerceptronNode hiddenNodes[16],
		PerceptronNode outputNodes[48])
{
	self->numInputNodes = numInputNodes;
	self->hiddenNodes = hiddenNodes;
	self->outputNodes = outputNodes;
}
 
void Perceptron_init(Perceptron *self, int numInputNodes)
{
	int i = 0;
	PerceptronNode_init(PerceptronNode *self, int numInputNodes);
 
	self->numInputNodes = numInputNodes;
	for(i = 0; i < 16; i++)
		self->hiddenNodes[i] = hiddenNodes[i];
	for(i = 0; i < 48; i++)
		self->outputNodes[i] = outputNodes[i];
}
 
void Perceptron_init(Perceptron *self, int numInputNodes, ranluxcl_state_t *rst)
{
}
 
/**
 * O(numHiddenNodes + numOutputNodes)
 */
void Perceptron_calc(Perceptron *self, float16 x, int outNodeIndex, float y[3])
{
	int i = 0;
	float16 hiddenNodeVals;
 
	// Calculate hidden nodes followed by output nodes
	hiddenNodeVals.s0 = PerceptronNode_calc(self->hiddenNodes[0], x);
	hiddenNodeVals.s1 = PerceptronNode_calc(self->hiddenNodes[1], x);
	hiddenNodeVals.s2 = PerceptronNode_calc(self->hiddenNodes[2], x);
	hiddenNodeVals.s3 = PerceptronNode_calc(self->hiddenNodes[3], x);
	hiddenNodeVals.s4 = PerceptronNode_calc(self->hiddenNodes[4], x);
	hiddenNodeVals.s5 = PerceptronNode_calc(self->hiddenNodes[5], x);
	hiddenNodeVals.s6 = PerceptronNode_calc(self->hiddenNodes[6], x);
	hiddenNodeVals.s7 = PerceptronNode_calc(self->hiddenNodes[7], x);
	hiddenNodeVals.s8 = PerceptronNode_calc(self->hiddenNodes[8], x);
	hiddenNodeVals.s9 = PerceptronNode_calc(self->hiddenNodes[9], x);
	hiddenNodeVals.sA = PerceptronNode_calc(self->hiddenNodes[10], x);
	hiddenNodeVals.sB = PerceptronNode_calc(self->hiddenNodes[11], x);
	hiddenNodeVals.sC = PerceptronNode_calc(self->hiddenNodes[12], x);
	hiddenNodeVals.sD = PerceptronNode_calc(self->hiddenNodes[13], x);
	hiddenNodeVals.sE = PerceptronNode_calc(self->hiddenNodes[14], x);
	hiddenNodeVals.sF = PerceptronNode_calc(self->hiddenNodes[15], x);
	y[0] = PerceptronNode_calc(self->outputNodes[outNodeIndex], hiddenNodeVals);
	y[1] = PerceptronNode_calc(self->outputNodes[outNodeIndex + 1],
			hiddenNodeVals);
	y[2] = PerceptronNode_calc(self->outputNodes[outNodeIndex + 2],
			hiddenNodeVals);
}
 
/**
 * O(numHiddenNodes + numOutputNodes)
 */
void Perceptron_calc(Perceptron *self, float16 x, float y[48])
{
	int i = 0;
	float16 hiddenNodeVals;
 
	// Calculate hidden nodes followed by output nodes
	hiddenNodeVals.s0 = PerceptronNode_calc(self->hiddenNodes[0], x);
	hiddenNodeVals.s1 = PerceptronNode_calc(self->hiddenNodes[1], x);
	hiddenNodeVals.s2 = PerceptronNode_calc(self->hiddenNodes[2], x);
	hiddenNodeVals.s3 = PerceptronNode_calc(self->hiddenNodes[3], x);
	hiddenNodeVals.s4 = PerceptronNode_calc(self->hiddenNodes[4], x);
	hiddenNodeVals.s5 = PerceptronNode_calc(self->hiddenNodes[5], x);
	hiddenNodeVals.s6 = PerceptronNode_calc(self->hiddenNodes[6], x);
	hiddenNodeVals.s7 = PerceptronNode_calc(self->hiddenNodes[7], x);
	hiddenNodeVals.s8 = PerceptronNode_calc(self->hiddenNodes[8], x);
	hiddenNodeVals.s9 = PerceptronNode_calc(self->hiddenNodes[9], x);
	hiddenNodeVals.sA = PerceptronNode_calc(self->hiddenNodes[10], x);
	hiddenNodeVals.sB = PerceptronNode_calc(self->hiddenNodes[11], x);
	hiddenNodeVals.sC = PerceptronNode_calc(self->hiddenNodes[12], x);
	hiddenNodeVals.sD = PerceptronNode_calc(self->hiddenNodes[13], x);
	hiddenNodeVals.sE = PerceptronNode_calc(self->hiddenNodes[14], x);
	hiddenNodeVals.sF = PerceptronNode_calc(self->hiddenNodes[15], x);
	for(i = 0; i < 48; i++)
		y[i] = PerceptronNode_calc(self->outputNodes[i], hiddenNodeVals);
}

GMM.cl

/*
 * GMM.cl
 *
 *  Created on: Jan 7, 2012
 *      Author: scannon
 */
 
#include "ranluxcl.cl"
#include "Perceptron.cl"
 
 
typedef struct __GMM__
{
	// 1 input per conditional variable
	// 3 outputs per normal:  mean, standard deviation, and weight
	Perceptron *perceptron;
	__GMM__ *conditional;
	int numConditionals;
	char init;
	float value;
} GMM;
 
 
char GMM__init__(GMM self, Perceptron *perceptron, GMM *conditionals, int numConditionals)
{
}
 
/**
 * O(numOutputNodes*O(normal())/3 + numHiddenNodes*(numInputs + numOutputNodes))
 */
float GMM_sample(GMM self, ranluxcl_state_t *rst, float4 *sample)
{
	float normParams[48];
	int outNodeIndex = 0;
	float4 ranluxRet = (0, 0, 0, 0);
	float mean = 0.0f;
	float sd = 0.0f;
	float weight = 0.0f;
	float weightSum = 0.0f;
	int i = 0;
 
	// Determine weights, means, and standard deviations from conditionals
	Perceptron_calc(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];
	*sample = ranluxcl32norm(rst)*sd + mean;
}

test.cl

/*
 * test.cl
 *
 *  Created on: Jan 22, 2012
 *      Author: scannon
 */
 
#ifndef TEST_CL
#define TEST_CL
 
 
 char test_PerceptronNode()
{
	//PerceptronNode pNode1 = PerceptronNode();
	PerceptronNode pNode1;
	float16 w1 = (float16)(0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f);
	float16 i1 = (float16)(1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f);
	PerceptronNode_init(&pNode1, w1);
	if(PerceptronNode_calc(&pNode1, i1) != 0.5f)
		return 1;
 
	//PerceptronNode pNode2 = PerceptronNode();
	PerceptronNode pNode2;
	float16 w2 = (float16)(1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 1.0f);
	float16 i2 = (float16)(0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f);
	PerceptronNode_init(&pNode2, w2);
	if(PerceptronNode_calc(&pNode2, i2) != 0.5f)
		return 2;
 
	//PerceptronNode pNode3 = PerceptronNode();
	PerceptronNode pNode3;
	float16 w3 = (float16)(1.0f, 1.0f, 1.0f, 1.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f);
	float16 i3 = (float16)(1.0f, 1.0f, 1.0f, 1.0f, 1.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f);
	PerceptronNode_init(&pNode3, w3);
	if(PerceptronNode_calc(&pNode3, i3) - 0.982f > 0.001f ||
	   PerceptronNode_calc(&pNode3, i3) - 0.982f < -0.001f)
		return 3;
 
	//PerceptronNode pNode4 = PerceptronNode();
	PerceptronNode pNode4;
	float16 w4 = (float16)(1.0f, 3.0f, 2.0f, 4.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f);
	float16 i4 = (float16)(1.0f, 3.0f, 5.0f, 7.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f);
	PerceptronNode_init(&pNode4, w4);
	if(PerceptronNode_calc(&pNode4, i4) - 0.999f > 0.001f ||
	   PerceptronNode_calc(&pNode4, i4) - 0.999f < -0.001f)
		return 4;
 
	return 0;
}
 
//char test_Perceptron()
//{
//	Perceptron p1 = Perdecptron();
//	Perceptron_init(&p1);
//}
 
kernel void Kernel_Ranluxcl_Init(private uint ins,
                                 global ranluxcl_state_t *ranluxcltab)
{
	ranluxcl_initialization(ins, ranluxcltab);
}
 
__kernel void test_kernel(__global const float *a,
                          __global const float *b,
                          private uint ins,
                          global ranluxcl_state_t *ranluxcltab,
                          __global char *result)
{
	int gid = get_global_id(0);
	ranluxcl_state_t ranluxclstate;
 
	//Download state
	ranluxcl_download_seed(&ranluxclstate, ranluxcltab);
 
	//Generate a float4 with each component on (0,1),
	//end points not included. We can call ranluxcl as many
	//times as we like until we upload the state again.
	float4 randomnr = ranluxcl32(&ranluxclstate);
 
	result[gid] = 0;
	result[gid] = test_PerceptronNode( );
	if(!result[gid])
		return;
 
	//Upload state again
	ranluxcl_upload_seed(&ranluxclstate, ranluxcltab);
 
}
 
#endif	// TEST_CL