andrew@2: /*
andrew@2:  *  bayesianArray.cpp
andrew@2:  *  bayesianTest5
andrew@2:  *
andrew@2:  *  Created by Andrew Robertson on 08/05/2010.
andrew@2:  *  Copyright 2010 __MyCompanyName__. All rights reserved.
andrew@2:  *
andrew@2:  */
andrew@2: 
andrew@2: #include "bayesianArray.h"
andrew@2: #include "math.h"
andrew@2: #include "ofMain.h"
andrew@2: 
andrew@2: bayesianArray::bayesianArray(){
andrew@2: 	likelihoodNoise = 0.5;
andrew@6: 	likelihoodMean = arraySize/2;
andrew@6: 	likelihoodStdDev = arraySize / 12;
andrew@2: initialiseArray();
andrew@2: }
andrew@2: 
andrew@2: void bayesianArray::initialiseArray(){
andrew@2: 
andrew@2: 	//maximumIndex = 12;//change this
andrew@6: 	setGaussianPrior(arraySize/2, arraySize/1);	
andrew@6: 	setGaussianLikelihood(arraySize/2, arraySize/1);//likelihoodMean, likelihoodStdDev);
andrew@2: 	
andrew@2: 	calculatePosterior();
andrew@2: 	renormalisePosterior();
andrew@2: 	posteriorDecayRate = 0.06;
andrew@2: 	
andrew@2: 	eighthNoteProportion = 0.35;//must be less than 0.5 to discriminate - was 0.4
andrew@2: 	earlySixteenthNoteProportion = 0;
andrew@2: 	lateSixteenthNoteProportion = 0;
andrew@2: 	decayNoiseAmount = 0.1;
andrew@6: 	decayNoiseStdDev = arraySize/24;
andrew@2: 	standardDeviation = likelihoodStdDev;
andrew@6: 	setDecayNoiseGaussian(arraySize/2, decayNoiseStdDev);
andrew@2: 	
andrew@2: 	setGaussianLikelihood(likelihoodMean, likelihoodStdDev);
andrew@2: }
andrew@2: 
andrew@2: 
andrew@2: void bayesianArray::setGaussianPrior(float mean, float StdDev){
andrew@2: int i;
andrew@6: 	for (i=0;i<arraySize;i++){
andrew@2: 	prior[i] = (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev));
andrew@2: 	//posterior[i] = prior[i];
andrew@2: 	}
andrew@2: }
andrew@2: 
andrew@2: void bayesianArray::setGaussianPosterior(float mean, float StdDev){
andrew@2: int i;
andrew@6: 	for (i=0;i<arraySize;i++){
andrew@2: 	posterior[i] = (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev));
andrew@2: 	}
andrew@2: }
andrew@2: 
andrew@2: void bayesianArray::setGaussianLikelihoodForBeats(float mean, float StdDev){
andrew@2: 	//this has eighth and sixteenth positions included
andrew@2: 	
andrew@6: if (mean >= 0 && mean <= arraySize){
andrew@2: int i;	float eighthDifference;
andrew@4: 	//main beat position plus these three:
andrew@6: int eighthPosition = ((int)mean + arraySize/2)%arraySize;
andrew@6: int earlySixteenthPosition = ((int)mean + (3*arraySize/4))%arraySize;;
andrew@6: int lateSixteenthPosition = ((int)mean + (arraySize/4))%arraySize;;
andrew@2: 
andrew@2: float mainDifference, sixteenthDifference;
andrew@2: float gaussianProportion = 1 - likelihoodNoise;
andrew@2: float mainProportion = (1 - eighthNoteProportion - earlySixteenthNoteProportion - lateSixteenthNoteProportion);
andrew@2: 
andrew@6: 	for (i=0;i < arraySize;i++){
andrew@2: 	
andrew@6: 	mainDifference = min( fabs(i-mean) , (double)(i + arraySize - mean));
andrew@2: 	likelihood[i] = gaussianProportion * mainProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(mainDifference)*(mainDifference)/(2*StdDev*StdDev)) ;
andrew@2: 	
andrew@6: 	eighthDifference = min( abs(i - eighthPosition) , i + arraySize - eighthPosition);
andrew@6: 	eighthDifference = min(eighthDifference , (float)(arraySize + eighthPosition - i ));
andrew@2: 	//for e.g. +0.43, or -0.47 we require the gaussian around the half note too
andrew@2: 	likelihood[i] += gaussianProportion * eighthNoteProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(eighthDifference)*(eighthDifference)/(2*StdDev*StdDev)) ;
andrew@2: 	
andrew@6: 	sixteenthDifference = min( abs(i - earlySixteenthPosition) , i + arraySize - earlySixteenthPosition);
andrew@6: 	sixteenthDifference = min(sixteenthDifference , (float)(arraySize + earlySixteenthPosition - i ));
andrew@2: 	//for e.g. +0.43, or -0.47 we require the gaussian around the half note too
andrew@2: 	likelihood[i] += gaussianProportion * earlySixteenthNoteProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(sixteenthDifference)*(sixteenthDifference)/(2*StdDev*StdDev)) ;
andrew@2: 	
andrew@6: 	sixteenthDifference = min( abs(i - lateSixteenthPosition) , i + arraySize - lateSixteenthPosition);
andrew@6: 	sixteenthDifference = min(sixteenthDifference , (float)(arraySize + lateSixteenthPosition - i ));
andrew@2: 	//for e.g. +0.43, or -0.47 we require the gaussian around the half note too
andrew@2: 	likelihood[i] += gaussianProportion * lateSixteenthNoteProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(sixteenthDifference)*(sixteenthDifference)/(2*StdDev*StdDev)) ;
andrew@2: 	
andrew@2: 
andrew@2: 
andrew@6: 	likelihood[i] += (likelihoodNoise / arraySize);
andrew@2: 	//likelihood[i] = (float) max(gaussianProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev)) , 
andrew@6: 	//(double) (likelihoodNoise / arraySize) );
andrew@2: 	}
andrew@6: //	renormaliseArray(&likelihood[0], arraySize);
andrew@2: }//end if mean within limits
andrew@2: }
andrew@2: 
andrew@2: 
andrew@2: void bayesianArray::setGaussianLikelihood(float mean, float StdDev){
andrew@6: if (mean >= 0 && mean <= arraySize){
andrew@2: int i;	float eighthDifference;
andrew@6: int eighthPosition = ((int)mean + arraySize/2)%arraySize;
andrew@2: float mainDifference;
andrew@2: float gaussianProportion = 1 - likelihoodNoise;
andrew@2: 
andrew@6: 	for (i=0;i < arraySize;i++){
andrew@2: 	
andrew@6: 	mainDifference = min( fabs(i-mean) , (double)(i + arraySize - mean));
andrew@2: 	//without * (1 - eighthNoteProportion)
andrew@2: 	likelihood[i] = gaussianProportion  * (1/(StdDev*sqrt(2*PI)))*exp(-1*(mainDifference)*(mainDifference)/(2*StdDev*StdDev)) ;
andrew@2: 		
andrew@6: 	likelihood[i] += (likelihoodNoise / arraySize);
andrew@2: 	//likelihood[i] = (float) max(gaussianProportion * (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev)) , 
andrew@6: 	//(double) (likelihoodNoise / arraySize) );
andrew@2: 	}
andrew@6: //	renormaliseArray(&likelihood[0], arraySize);
andrew@2: }//end if mean within limits
andrew@2: }
andrew@2: 
andrew@2: void bayesianArray::calculatePosterior(){
andrew@2: int i;
andrew@6: 	for (i=0;i < arraySize;i++){
andrew@2: 	posterior[i] = likelihood[i] * prior[i];
andrew@2: 	}
andrew@2: 	//renormalisePosterior();
andrew@2: }
andrew@2: 
andrew@2: 
andrew@2: float bayesianArray::getMaximum(float *ptr, int length){
andrew@2: int i;
andrew@2: float max = 0;
andrew@2: 	for (i=0;i < length;i++){	
andrew@2: 		if (*(ptr+i)>max)
andrew@2: 		max = *(ptr+i);
andrew@2: 	}
andrew@2: 	maximumValue = max;
andrew@2: 	return max;
andrew@2: }
andrew@2: 
andrew@2: float* bayesianArray::getMaximumEstimate(float *ptr, int length){
andrew@2: float returnArray[2];
andrew@2: int i;
andrew@2: float max = 0;
andrew@2: maximumIndex = 0;
andrew@2: 	for (i=0;i < length;i++){	
andrew@2: 		if (*(ptr+i)>max){
andrew@2: 		max = *(ptr+i);
andrew@2: 		maximumIndex = i;
andrew@2: 		}
andrew@2: 	}
andrew@2: 	returnArray[0] = max;
andrew@2: 	returnArray[1] = maximumIndex;
andrew@2: 	maximumValue = max;
andrew@2: 	return &returnArray[0];
andrew@2: }
andrew@2: 
andrew@2: 
andrew@2: 
andrew@2: double bayesianArray::getIntegratedEstimateIndex(){
andrew@2: 	int i;
andrew@2: 	float integratedQuantity = 0;
andrew@2: 	float integratedTotal = 0;
andrew@2: 	double integratedIndex = 0;
andrew@6: 	for (i=0;i < arraySize;i++){	
andrew@2: 		integratedQuantity += posterior[i];//the values of the probability distribution
andrew@2: 		integratedTotal += i*posterior[i];
andrew@2: 	}
andrew@2: 	if (integratedQuantity > 0){
andrew@2: 		integratedIndex = integratedTotal / integratedQuantity;
andrew@2: 	}
andrew@2: 	integratedEstimate = (float) integratedIndex;
andrew@2: 	return integratedIndex;
andrew@2: }
andrew@2: 
andrew@2: 
andrew@2: double bayesianArray::calculateStandardDeviation(){
andrew@2: 
andrew@2: 	double total = 0;
andrew@2: 	double pdfSum;
andrew@2: 	double variance = 0;
andrew@6: 	for (int i=0;i < arraySize;i++){	
andrew@2: //*posterior[i] * 		
andrew@2: 		total += posterior[i] * (i - integratedEstimate) * (i - integratedEstimate);//the values of the probability distribution
andrew@2: 		pdfSum += posterior[i];
andrew@2: 	}
andrew@2: 	
andrew@2: 	if (pdfSum > 0)
andrew@2: 	variance = total / pdfSum;
andrew@2: 	else 
andrew@6: 	variance = arraySize;
andrew@2: 
andrew@2: 	standardDeviation = sqrt(variance);
andrew@2: 	return standardDeviation;
andrew@2: }
andrew@2: 
andrew@2: 
andrew@2: 
andrew@2: void bayesianArray::renormaliseArray(float *ptr, int length){
andrew@2: int i;
andrew@2: float totalArea = 0;
andrew@2: 	for (i=0;i < length;i++){
andrew@2: 		totalArea += *(ptr+i);
andrew@2: 	}
andrew@2: 	
andrew@2: 	for (i=0;i < length;i++){
andrew@2: 		*(ptr+i) /= totalArea;
andrew@2: 	}
andrew@2: 	
andrew@2: }
andrew@2: 
andrew@2: void bayesianArray::resetPrior(){
andrew@2: 	int i;
andrew@6: 	for (i=0;i<arraySize;i++){
andrew@2: 	prior[i] = posterior[i];
andrew@2: 	}
andrew@2: }
andrew@2: 
andrew@2: void bayesianArray::renormalisePrior(){
andrew@2: 	int i;
andrew@2: 	float totalArea = 0;
andrew@6: 	for (i=0;i < arraySize;i++){
andrew@2: 		totalArea += prior[i];
andrew@2: 	}
andrew@2: 	
andrew@6: 	for (i=0;i < arraySize;i++){
andrew@2: 		prior[i] /= totalArea;
andrew@2: 	}
andrew@2: }
andrew@2: 
andrew@2: void bayesianArray::renormalisePosterior(){
andrew@2: int i;
andrew@2: float totalArea = 0;
andrew@6: 	for (i=0;i < arraySize;i++){
andrew@2: 		totalArea += posterior[i];
andrew@2: 	}
andrew@2: 	
andrew@6: 	for (i=0;i < arraySize;i++){
andrew@2: 		posterior[i] /= totalArea;
andrew@2: 	}
andrew@2: }
andrew@2: 
andrew@2: void bayesianArray::decayPosterior(){
andrew@2: 	float *pointer;
andrew@6: 	pointer = getMaximumEstimate(&posterior[0], arraySize);	
andrew@2: 	 float maximum;
andrew@2: 	 maximum = *pointer;
andrew@2: 	int i;
andrew@6: 	for (i=0;i<arraySize;i++){
andrew@2: 	posterior[i] += (maximum - posterior[i]) * posteriorDecayRate * 0.01;;//usded to be * maximum not minus value
andrew@2: 	}
andrew@2: 	maximumIndex = *(pointer+1);
andrew@2: }
andrew@2: 
andrew@2: void bayesianArray::setDecayNoiseGaussian(float mean, float StdDev){
andrew@2: int i;
andrew@6: 	for (i=0;i<arraySize;i++){
andrew@2: 	decayNoiseArray[i] = (1/(StdDev*sqrt(2*PI)))*exp(-1*(i-mean)*(i-mean)/(2*StdDev*StdDev));
andrew@2: 	}
andrew@2: }
andrew@2: 
andrew@2: void bayesianArray::decayPosteriorWithGaussianNoise(){
andrew@2: 
andrew@2: 	int i;
andrew@6: 	float currentMaximum = getMaximum(&posterior[0], arraySize);
andrew@6: 	for (i=0;i<arraySize;i++){
andrew@6: 	posterior[i] += decayNoiseArray[(i - (int)maximumIndex + ((3*arraySize)/2)) % arraySize] * currentMaximum * decayNoiseAmount;
andrew@2: 	//posteriorDecayRate * 0.01;;//usded to be * maximum not minus value
andrew@2: 	}
andrew@2: 
andrew@2: }
andrew@2: 
andrew@2: void bayesianArray::resetMaximumPosterior(){
andrew@2: int i;
andrew@2: float max = 0;
andrew@6: 	for (i=0;i < arraySize;i++){	
andrew@2: 		if (posterior[i]>max){
andrew@2: 		maximumIndex = i;
andrew@2: 		max = posterior[i];
andrew@2: 		}
andrew@2: 	}
andrew@2: }
andrew@2: 
andrew@2: void bayesianArray::translateDistribution(int translationIndex){
andrew@2: int tmpIndex;
andrew@2: 		//copy array
andrew@2: 		int i;
andrew@6: 		for (i=0;i < arraySize;i++){
andrew@2: 		tempPosteriorArray[i] =	posterior[i] ;
andrew@2: 		}
andrew@2: 		//translate values
andrew@6: 	for (i=0;i < arraySize;i++){
andrew@6: 	tmpIndex = (i + translationIndex + arraySize)%arraySize;
andrew@2: 	posterior[tmpIndex] = tempPosteriorArray[i]; 
andrew@2: 	}
andrew@2: 	//now delete tmp array
andrew@2: }
andrew@2: 
andrew@2: double bayesianArray::getKLdivergence(){
andrew@2: 	double KLsum = 0;
andrew@2: 	//take no chances - renormalise both prior and posterior
andrew@2: 	renormalisePosterior();
andrew@2: 	renormalisePrior();
andrew@6: 	for (int i = 0;i < arraySize;i++){
andrew@2: 		if (posterior[i] > 0 && prior[i] > 0){
andrew@2: 			KLsum += (posterior[i]*log(posterior[i]/prior[i]));
andrew@2: 		}
andrew@2: 	}
andrew@2: 	return KLsum;
andrew@2: }
andrew@2: 
andrew@3: double bayesianArray::getEntropyOfPosterior(){
andrew@3: 	double entropy = 0;
andrew@4: 	//make sure normalised? (it is)
andrew@6: 	for (int i = 0;i < arraySize;i++){
andrew@3: 		if (posterior[i] > 0){
andrew@3: 			entropy += posterior[i]*log(posterior[i]);
andrew@3: 		}
andrew@3: 	}
andrew@3: 	return -1*entropy;
andrew@3: }
andrew@3: 
andrew@11: double bayesianArray::getEntropyOfPrior(){
andrew@11: 	double entropy = 0;
andrew@11: 	//make sure normalised? (it is)
andrew@11: 	for (int i = 0;i < arraySize;i++){
andrew@11: 		if (posterior[i] > 0){
andrew@11: 			entropy += prior[i]*log(prior[i]);
andrew@11: 		}
andrew@11: 	}
andrew@11: 	return -1*entropy;
andrew@11: }
andrew@3: 
andrew@11: