view newOFsrc/bayesianArray.cpp @ 6:fc095148e2a8

sending out the KL information that cognitive dynamic paper will be using

/*
 *  bayesianArray.cpp
 *  bayesianTest5
 *
 *  Created by Andrew Robertson on 08/05/2010.
 *  Copyright 2010 __MyCompanyName__. All rights reserved.
 *
 */

#include "bayesianArray.h"
#include "math.h"
#include "ofMain.h"

bayesianArray::bayesianArray(){
	likelihoodNoise = 0.5;
	likelihoodMean = arraySize/2;
	likelihoodStdDev = arraySize / 12;
initialiseArray();
}

void bayesianArray::initialiseArray(){

	//maximumIndex = 12;//change this
	setGaussianPrior(arraySize/2, arraySize/1);	
	setGaussianLikelihood(arraySize/2, arraySize/1);//likelihoodMean, likelihoodStdDev);
	
	calculatePosterior();
	renormalisePosterior();
	posteriorDecayRate = 0.06;
	
	eighthNoteProportion = 0.35;//must be less than 0.5 to discriminate - was 0.4
	earlySixteenthNoteProportion = 0;
	lateSixteenthNoteProportion = 0;
	decayNoiseAmount = 0.1;
	decayNoiseStdDev = arraySize/24;
	standardDeviation = likelihoodStdDev;
	setDecayNoiseGaussian(arraySize/2, decayNoiseStdDev);
	
	setGaussianLikelihood(likelihoodMean, likelihoodStdDev);
}


void bayesianArray::setGaussianPrior(float mean, float StdDev){
int i;
	for (i=0;i<arraySize;i++){
	prior[i] = (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev));
	//posterior[i] = prior[i];
	}
}

void bayesianArray::setGaussianPosterior(float mean, float StdDev){
int i;
	for (i=0;i<arraySize;i++){
	posterior[i] = (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev));
	}
}

void bayesianArray::setGaussianLikelihoodForBeats(float mean, float StdDev){
	//this has eighth and sixteenth positions included
	
if (mean >= 0 && mean <= arraySize){
int i;	float eighthDifference;
	//main beat position plus these three:
int eighthPosition = ((int)mean + arraySize/2)%arraySize;
int earlySixteenthPosition = ((int)mean + (3*arraySize/4))%arraySize;;
int lateSixteenthPosition = ((int)mean + (arraySize/4))%arraySize;;

float mainDifference, sixteenthDifference;
float gaussianProportion = 1 - likelihoodNoise;
float mainProportion = (1 - eighthNoteProportion - earlySixteenthNoteProportion - lateSixteenthNoteProportion);

	for (i=0;i < arraySize;i++){
	
	mainDifference = min( fabs(i-mean) , (double)(i + arraySize - mean));
	likelihood[i] = gaussianProportion * mainProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(mainDifference)*(mainDifference)/(2*StdDev*StdDev)) ;
	
	eighthDifference = min( abs(i - eighthPosition) , i + arraySize - eighthPosition);
	eighthDifference = min(eighthDifference , (float)(arraySize + eighthPosition - i ));
	//for e.g. +0.43, or -0.47 we require the gaussian around the half note too
	likelihood[i] += gaussianProportion * eighthNoteProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(eighthDifference)*(eighthDifference)/(2*StdDev*StdDev)) ;
	
	sixteenthDifference = min( abs(i - earlySixteenthPosition) , i + arraySize - earlySixteenthPosition);
	sixteenthDifference = min(sixteenthDifference , (float)(arraySize + earlySixteenthPosition - i ));
	//for e.g. +0.43, or -0.47 we require the gaussian around the half note too
	likelihood[i] += gaussianProportion * earlySixteenthNoteProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(sixteenthDifference)*(sixteenthDifference)/(2*StdDev*StdDev)) ;
	
	sixteenthDifference = min( abs(i - lateSixteenthPosition) , i + arraySize - lateSixteenthPosition);
	sixteenthDifference = min(sixteenthDifference , (float)(arraySize + lateSixteenthPosition - i ));
	//for e.g. +0.43, or -0.47 we require the gaussian around the half note too
	likelihood[i] += gaussianProportion * lateSixteenthNoteProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(sixteenthDifference)*(sixteenthDifference)/(2*StdDev*StdDev)) ;
	


	likelihood[i] += (likelihoodNoise / arraySize);
	//likelihood[i] = (float) max(gaussianProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev)) , 
	//(double) (likelihoodNoise / arraySize) );
	}
//	renormaliseArray(&likelihood[0], arraySize);
}//end if mean within limits
}


void bayesianArray::setGaussianLikelihood(float mean, float StdDev){
if (mean >= 0 && mean <= arraySize){
int i;	float eighthDifference;
int eighthPosition = ((int)mean + arraySize/2)%arraySize;
float mainDifference;
float gaussianProportion = 1 - likelihoodNoise;

	for (i=0;i < arraySize;i++){
	
	mainDifference = min( fabs(i-mean) , (double)(i + arraySize - mean));
	//without * (1 - eighthNoteProportion)
	likelihood[i] = gaussianProportion  * (1/(StdDev*sqrt(2*PI)))*exp(-1*(mainDifference)*(mainDifference)/(2*StdDev*StdDev)) ;
		
	likelihood[i] += (likelihoodNoise / arraySize);
	//likelihood[i] = (float) max(gaussianProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev)) , 
	//(double) (likelihoodNoise / arraySize) );
	}
//	renormaliseArray(&likelihood[0], arraySize);
}//end if mean within limits
}

void bayesianArray::calculatePosterior(){
int i;
	for (i=0;i < arraySize;i++){
	posterior[i] = likelihood[i] * prior[i];
	}
	//renormalisePosterior();
}


float bayesianArray::getMaximum(float *ptr, int length){
int i;
float max = 0;
	for (i=0;i < length;i++){	
		if (*(ptr+i)>max)
		max = *(ptr+i);
	}
	maximumValue = max;
	return max;
}

float* bayesianArray::getMaximumEstimate(float *ptr, int length){
float returnArray[2];
int i;
float max = 0;
maximumIndex = 0;
	for (i=0;i < length;i++){	
		if (*(ptr+i)>max){
		max = *(ptr+i);
		maximumIndex = i;
		}
	}
	returnArray[0] = max;
	returnArray[1] = maximumIndex;
	maximumValue = max;
	return &returnArray[0];
}



double bayesianArray::getIntegratedEstimateIndex(){
	int i;
	float integratedQuantity = 0;
	float integratedTotal = 0;
	double integratedIndex = 0;
	for (i=0;i < arraySize;i++){	
		integratedQuantity += posterior[i];//the values of the probability distribution
		integratedTotal += i*posterior[i];
	}
	if (integratedQuantity > 0){
		integratedIndex = integratedTotal / integratedQuantity;
	}
	integratedEstimate = (float) integratedIndex;
	return integratedIndex;
}


double bayesianArray::calculateStandardDeviation(){

	double total = 0;
	double pdfSum;
	double variance = 0;
	for (int i=0;i < arraySize;i++){	
//*posterior[i] * 		
		total += posterior[i] * (i - integratedEstimate) * (i - integratedEstimate);//the values of the probability distribution
		pdfSum += posterior[i];
	}
	
	if (pdfSum > 0)
	variance = total / pdfSum;
	else 
	variance = arraySize;

	standardDeviation = sqrt(variance);
	return standardDeviation;
}



void bayesianArray::renormaliseArray(float *ptr, int length){
int i;
float totalArea = 0;
	for (i=0;i < length;i++){
		totalArea += *(ptr+i);
	}
	
	for (i=0;i < length;i++){
		*(ptr+i) /= totalArea;
	}
	
}

void bayesianArray::resetPrior(){
	int i;
	for (i=0;i<arraySize;i++){
	prior[i] = posterior[i];
	}
}

void bayesianArray::renormalisePrior(){
	int i;
	float totalArea = 0;
	for (i=0;i < arraySize;i++){
		totalArea += prior[i];
	}
	
	for (i=0;i < arraySize;i++){
		prior[i] /= totalArea;
	}
}

void bayesianArray::renormalisePosterior(){
int i;
float totalArea = 0;
	for (i=0;i < arraySize;i++){
		totalArea += posterior[i];
	}
	
	for (i=0;i < arraySize;i++){
		posterior[i] /= totalArea;
	}
}

void bayesianArray::decayPosterior(){
	float *pointer;
	pointer = getMaximumEstimate(&posterior[0], arraySize);	
	 float maximum;
	 maximum = *pointer;
	int i;
	for (i=0;i<arraySize;i++){
	posterior[i] += (maximum - posterior[i]) * posteriorDecayRate * 0.01;;//usded to be * maximum not minus value
	}
	maximumIndex = *(pointer+1);
}

void bayesianArray::setDecayNoiseGaussian(float mean, float StdDev){
int i;
	for (i=0;i<arraySize;i++){
	decayNoiseArray[i] = (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev));
	}
}

void bayesianArray::decayPosteriorWithGaussianNoise(){

	int i;
	float currentMaximum = getMaximum(&posterior[0], arraySize);
	for (i=0;i<arraySize;i++){
	posterior[i] += decayNoiseArray[(i - (int)maximumIndex + ((3*arraySize)/2)) % arraySize] * currentMaximum * decayNoiseAmount;
	//posteriorDecayRate * 0.01;;//usded to be * maximum not minus value
	}

}

void bayesianArray::resetMaximumPosterior(){
int i;
float max = 0;
	for (i=0;i < arraySize;i++){	
		if (posterior[i]>max){
		maximumIndex = i;
		max = posterior[i];
		}
	}
}

void bayesianArray::translateDistribution(int translationIndex){
int tmpIndex;
		//copy array
		int i;
		for (i=0;i < arraySize;i++){
		tempPosteriorArray[i] =	posterior[i] ;
		}
		//translate values
	for (i=0;i < arraySize;i++){
	tmpIndex = (i + translationIndex + arraySize)%arraySize;
	posterior[tmpIndex] = tempPosteriorArray[i]; 
	}
	//now delete tmp array
}

double bayesianArray::getKLdivergence(){
	double KLsum = 0;
	//take no chances - renormalise both prior and posterior
	renormalisePosterior();
	renormalisePrior();
	for (int i = 0;i < arraySize;i++){
		if (posterior[i] > 0 && prior[i] > 0){
			KLsum += (posterior[i]*log(posterior[i]/prior[i]));
		}
	}
	return KLsum;
}

double bayesianArray::getEntropyOfPosterior(){
	double entropy = 0;
	//make sure normalised? (it is)
	for (int i = 0;i < arraySize;i++){
		if (posterior[i] > 0){
			entropy += posterior[i]*log(posterior[i]);
		}
	}
	return -1*entropy;
}